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 "CappedDruckerPragerCosseratStressUpdate.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", CappedDruckerPragerCosseratStressUpdate);
      13             : 
      14             : InputParameters
      15         168 : CappedDruckerPragerCosseratStressUpdate::validParams()
      16             : {
      17         168 :   InputParameters params = CappedDruckerPragerStressUpdate::validParams();
      18         168 :   params.addClassDescription("Capped Drucker-Prager plasticity stress calculator for the Cosserat "
      19             :                              "situation where the host medium (ie, the limit where all Cosserat "
      20             :                              "effects are zero) is isotropic.  Note that the return-map flow rule "
      21             :                              "uses an isotropic elasticity tensor built with the 'host' properties "
      22             :                              "defined by the user.");
      23         336 :   params.addRequiredRangeCheckedParam<Real>("host_youngs_modulus",
      24             :                                             "host_youngs_modulus>0",
      25             :                                             "Young's modulus for the isotropic host medium");
      26         336 :   params.addRequiredRangeCheckedParam<Real>("host_poissons_ratio",
      27             :                                             "host_poissons_ratio>=0 & host_poissons_ratio<0.5",
      28             :                                             "Poisson's ratio for the isotropic host medium");
      29         168 :   return params;
      30           0 : }
      31             : 
      32         126 : CappedDruckerPragerCosseratStressUpdate::CappedDruckerPragerCosseratStressUpdate(
      33         126 :     const InputParameters & parameters)
      34             :   : CappedDruckerPragerStressUpdate(parameters),
      35         504 :     _shear(getParam<Real>("host_youngs_modulus") /
      36         378 :            (2.0 * (1.0 + getParam<Real>("host_poissons_ratio"))))
      37             : {
      38         252 :   const Real young = getParam<Real>("host_youngs_modulus");
      39         252 :   const Real poisson = getParam<Real>("host_poissons_ratio");
      40         126 :   const Real lambda = young * poisson / ((1.0 + poisson) * (1.0 - 2.0 * poisson));
      41         126 :   _Ehost.fillFromInputVector({lambda, _shear}, RankFourTensor::symmetric_isotropic);
      42         126 : }
      43             : 
      44             : void
      45        7200 : CappedDruckerPragerCosseratStressUpdate::setEppEqq(const RankFourTensor & /*Eijkl*/,
      46             :                                                    Real & Epp,
      47             :                                                    Real & Eqq) const
      48             : {
      49        7200 :   Epp = _Ehost.sum3x3();
      50        7200 :   Eqq = _shear;
      51        7200 : }
      52             : 
      53             : void
      54        7200 : CappedDruckerPragerCosseratStressUpdate::setStressAfterReturn(const RankTwoTensor & stress_trial,
      55             :                                                               Real p_ok,
      56             :                                                               Real q_ok,
      57             :                                                               Real /*gaE*/,
      58             :                                                               const std::vector<Real> & /*intnl*/,
      59             :                                                               const yieldAndFlow & /*smoothed_q*/,
      60             :                                                               const RankFourTensor & /*Eijkl*/,
      61             :                                                               RankTwoTensor & stress) const
      62             : {
      63             :   // symm_stress is the symmetric part of the stress tensor.
      64             :   // symm_stress = (s_ij+s_ji)/2 + de_ij tr(stress) / 3
      65             :   //             = q / q_trial * (s_ij^trial+s_ji^trial)/2 + de_ij p / 3
      66             :   //             = q / q_trial * (symm_stress_ij^trial - de_ij tr(stress^trial) / 3) + de_ij p / 3
      67        7200 :   const Real p_trial = stress_trial.trace();
      68       14400 :   RankTwoTensor symm_stress = RankTwoTensor(RankTwoTensor::initIdentity) / 3.0 *
      69        7200 :                               (p_ok - (_in_q_trial == 0.0 ? 0.0 : p_trial * q_ok / _in_q_trial));
      70        7200 :   if (_in_q_trial > 0)
      71        7200 :     symm_stress += q_ok / _in_q_trial * 0.5 * (stress_trial + stress_trial.transpose());
      72        7200 :   stress = symm_stress + 0.5 * (stress_trial - stress_trial.transpose());
      73        7200 : }
      74             : 
      75             : void
      76         160 : CappedDruckerPragerCosseratStressUpdate::consistentTangentOperator(
      77             :     const RankTwoTensor & /*stress_trial*/,
      78             :     Real /*p_trial*/,
      79             :     Real /*q_trial*/,
      80             :     const RankTwoTensor & stress,
      81             :     Real /*p*/,
      82             :     Real q,
      83             :     Real gaE,
      84             :     const yieldAndFlow & smoothed_q,
      85             :     const RankFourTensor & Eijkl,
      86             :     bool compute_full_tangent_operator,
      87             :     RankFourTensor & cto) const
      88             : {
      89         160 :   if (!compute_full_tangent_operator)
      90             :   {
      91           0 :     cto = Eijkl;
      92           0 :     return;
      93             :   }
      94             : 
      95         160 :   RankFourTensor EAijkl;
      96         640 :   for (unsigned i = 0; i < _tensor_dimensionality; ++i)
      97        1920 :     for (unsigned j = 0; j < _tensor_dimensionality; ++j)
      98        5760 :       for (unsigned k = 0; k < _tensor_dimensionality; ++k)
      99       17280 :         for (unsigned l = 0; l < _tensor_dimensionality; ++l)
     100             :         {
     101       12960 :           cto(i, j, k, l) = 0.5 * (Eijkl(i, j, k, l) + Eijkl(j, i, k, l));
     102       12960 :           EAijkl(i, j, k, l) = 0.5 * (Eijkl(i, j, k, l) - Eijkl(j, i, k, l));
     103             :         }
     104             : 
     105             :   const RankTwoTensor s_over_q =
     106         320 :       (q == 0.0 ? RankTwoTensor()
     107         640 :                 : (0.5 * (stress + stress.transpose()) -
     108         800 :                    stress.trace() * RankTwoTensor(RankTwoTensor::initIdentity) / 3.0) /
     109             :                       q);
     110         160 :   const RankTwoTensor E_s_over_q = Eijkl.innerProductTranspose(s_over_q); // not symmetric in kl
     111             :   const RankTwoTensor Ekl =
     112         160 :       RankTwoTensor(RankTwoTensor::initIdentity).initialContraction(Eijkl); // symmetric in kl
     113             : 
     114         640 :   for (unsigned i = 0; i < _tensor_dimensionality; ++i)
     115        1920 :     for (unsigned j = 0; j < _tensor_dimensionality; ++j)
     116        5760 :       for (unsigned k = 0; k < _tensor_dimensionality; ++k)
     117       17280 :         for (unsigned l = 0; l < _tensor_dimensionality; ++l)
     118             :         {
     119       21600 :           cto(i, j, k, l) -= (i == j) * (1.0 / 3.0) *
     120       12960 :                              (Ekl(k, l) * (1.0 - _dp_dpt) + 0.5 * E_s_over_q(k, l) * (-_dp_dqt));
     121       12960 :           cto(i, j, k, l) -=
     122       12960 :               s_over_q(i, j) * (Ekl(k, l) * (-_dq_dpt) + 0.5 * E_s_over_q(k, l) * (1.0 - _dq_dqt));
     123             :         }
     124             : 
     125         160 :   if (smoothed_q.dg[1] != 0.0)
     126             :   {
     127         192 :     const RankFourTensor Tijab = _Ehost * (gaE / _Epp) * smoothed_q.dg[1] * d2qdstress2(stress);
     128          96 :     RankFourTensor inv = RankFourTensor(RankFourTensor::initIdentitySymmetricFour) + Tijab;
     129             :     try
     130             :     {
     131          96 :       inv = inv.transposeMajor().invSymm();
     132             :     }
     133           0 :     catch (const MooseException & e)
     134             :     {
     135             :       // Cannot form the inverse, so probably at some degenerate place in stress space.
     136             :       // Just return with the "best estimate" of the cto.
     137           0 :       mooseWarning("CappedDruckerPragerCosseratStressUpdate: Cannot invert 1+T in consistent "
     138             :                    "tangent operator computation at quadpoint ",
     139           0 :                    _qp,
     140             :                    " of element ",
     141           0 :                    _current_elem->id());
     142             :       return;
     143           0 :     }
     144          96 :     cto = (cto.transposeMajor() * inv).transposeMajor();
     145             :   }
     146         160 :   cto += EAijkl;
     147             : }