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 "CappedMohrCoulombCosseratStressUpdate.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("SolidMechanicsApp", CappedMohrCoulombCosseratStressUpdate);
      14             : 
      15             : InputParameters
      16         204 : CappedMohrCoulombCosseratStressUpdate::validParams()
      17             : {
      18         204 :   InputParameters params = CappedMohrCoulombStressUpdate::validParams();
      19         204 :   params.addClassDescription("Capped Mohr-Coulomb plasticity stress calculator for the Cosserat "
      20             :                              "situation where the host medium (ie, the limit where all Cosserat "
      21             :                              "effects are zero) is isotropic.  Note that the return-map flow rule "
      22             :                              "uses an isotropic elasticity tensor built with the 'host' properties "
      23             :                              "defined by the user.");
      24         408 :   params.addRequiredRangeCheckedParam<Real>("host_youngs_modulus",
      25             :                                             "host_youngs_modulus>0",
      26             :                                             "Young's modulus for the isotropic host medium");
      27         408 :   params.addRequiredRangeCheckedParam<Real>("host_poissons_ratio",
      28             :                                             "host_poissons_ratio>=0 & host_poissons_ratio<0.5",
      29             :                                             "Poisson's ratio for the isotropic host medium");
      30         204 :   return params;
      31           0 : }
      32             : 
      33         153 : CappedMohrCoulombCosseratStressUpdate::CappedMohrCoulombCosseratStressUpdate(
      34         153 :     const InputParameters & parameters)
      35             :   : CappedMohrCoulombStressUpdate(parameters),
      36         153 :     _host_young(getParam<Real>("host_youngs_modulus")),
      37         306 :     _host_poisson(getParam<Real>("host_poissons_ratio")),
      38         153 :     _host_E0011(_host_young * _host_poisson / (1.0 + _host_poisson) / (1.0 - 2.0 * _host_poisson)),
      39         153 :     _host_E0000(_host_E0011 + _host_young / (1.0 + _host_poisson)),
      40         306 :     _eigvals_scratch(_tensor_dimensionality)
      41             : {
      42         153 : }
      43             : 
      44             : void
      45        5904 : CappedMohrCoulombCosseratStressUpdate::preReturnMapV(
      46             :     const std::vector<Real> & /*trial_stress_params*/,
      47             :     const RankTwoTensor & stress_trial,
      48             :     const std::vector<Real> & /*intnl_old*/,
      49             :     const std::vector<Real> & /*yf*/,
      50             :     const RankFourTensor & /*Eijkl*/)
      51             : {
      52             :   mooseAssert(
      53             :       _eigvals_scratch.size() == _tensor_dimensionality,
      54             :       "_eigvals_scratch incorrectly sized in CappedMohCoulombCosseratStressUpdate::preReturnMapV");
      55        5904 :   stress_trial.symmetricEigenvaluesEigenvectors(_eigvals_scratch, _eigvecs);
      56        5904 :   _poissons_ratio = _host_poisson;
      57        5904 : }
      58             : 
      59             : void
      60        5904 : CappedMohrCoulombCosseratStressUpdate::setEffectiveElasticity(const RankFourTensor & /*Eijkl*/)
      61             : {
      62        5904 :   _Eij[0][0] = _Eij[1][1] = _Eij[2][2] = _host_E0000;
      63        5904 :   _Eij[0][1] = _Eij[1][0] = _Eij[0][2] = _Eij[2][0] = _Eij[1][2] = _Eij[2][1] = _host_E0011;
      64        5904 :   _En = _Eij[2][2];
      65        5904 :   const Real denom = _Eij[0][0] * (_Eij[0][0] + _Eij[0][1]) - 2 * Utility::pow<2>(_Eij[0][1]);
      66       23616 :   for (unsigned a = 0; a < _num_sp; ++a)
      67             :   {
      68       17712 :     _Cij[a][a] = (_Eij[0][0] + _Eij[0][1]) / denom;
      69       35424 :     for (unsigned b = 0; b < a; ++b)
      70       17712 :       _Cij[a][b] = _Cij[b][a] = -_Eij[0][1] / denom;
      71             :   }
      72        5904 : }
      73             : 
      74             : void
      75        5904 : CappedMohrCoulombCosseratStressUpdate::setStressAfterReturnV(
      76             :     const RankTwoTensor & stress_trial,
      77             :     const std::vector<Real> & stress_params,
      78             :     Real /*gaE*/,
      79             :     const std::vector<Real> & /*intnl*/,
      80             :     const yieldAndFlow & /*smoothed_q*/,
      81             :     const RankFourTensor & /*Eijkl*/,
      82             :     RankTwoTensor & stress) const
      83             : {
      84             :   // form the diagonal stress
      85        5904 :   stress = RankTwoTensor(stress_params[0], stress_params[1], stress_params[2], 0.0, 0.0, 0.0);
      86             :   // rotate to the original frame, to give the symmetric part of the stress
      87       11808 :   stress = _eigvecs * stress * (_eigvecs.transpose());
      88             :   // add the non-symmetric parts
      89        5904 :   stress += 0.5 * (stress_trial - stress_trial.transpose());
      90        5904 : }
      91             : 
      92             : void
      93         128 : CappedMohrCoulombCosseratStressUpdate::consistentTangentOperatorV(
      94             :     const RankTwoTensor & stress_trial,
      95             :     const std::vector<Real> & trial_stress_params,
      96             :     const RankTwoTensor & stress,
      97             :     const std::vector<Real> & stress_params,
      98             :     Real gaE,
      99             :     const yieldAndFlow & smoothed_q,
     100             :     const RankFourTensor & elasticity_tensor,
     101             :     bool compute_full_tangent_operator,
     102             :     const std::vector<std::vector<Real>> & dvar_dtrial,
     103             :     RankFourTensor & cto)
     104             : {
     105         128 :   CappedMohrCoulombStressUpdate::consistentTangentOperatorV(stress_trial,
     106             :                                                             trial_stress_params,
     107             :                                                             stress,
     108             :                                                             stress_params,
     109             :                                                             gaE,
     110             :                                                             smoothed_q,
     111             :                                                             elasticity_tensor,
     112             :                                                             compute_full_tangent_operator,
     113             :                                                             dvar_dtrial,
     114             :                                                             cto);
     115             : 
     116         128 :   if (!compute_full_tangent_operator)
     117           0 :     return;
     118             : 
     119             :   /**
     120             :    * Add the correction for the antisymmetric part of the elasticity
     121             :    * tensor.
     122             :    * CappedMohrCoulombStressUpdate computes
     123             :    * cto(i, j, k, l) = dstress(i, j)/dstrain(k, l)
     124             :    * and during the computations it explicitly performs certain
     125             :    * contractions that result in symmetry between i and j, and k and l,
     126             :    * viz, cto(i, j, k, l) = cto(j, i, k, l) = cto(i, j, l, k)
     127             :    * That is correct because that plasticity model is only valid for
     128             :    * symmetric stresses and strains.
     129             :    * CappedMohrCoulombCosseratStressUpdate does not include contributions from the
     130             :    * antisymmetric parts of stress (or strain), so the antisymmetric
     131             :    * parts of cto are just the antisymmetric parts of the elasticity
     132             :    * tensor, which must now get added to the cto computed by
     133             :    * CappedMohrCoulombStressUpdate
     134             :    */
     135         128 :   RankFourTensor anti;
     136         512 :   for (unsigned i = 0; i < _tensor_dimensionality; ++i)
     137        1536 :     for (unsigned j = 0; j < _tensor_dimensionality; ++j)
     138        4608 :       for (unsigned k = 0; k < _tensor_dimensionality; ++k)
     139       13824 :         for (unsigned l = 0; l < _tensor_dimensionality; ++l)
     140       10368 :           anti(i, j, k, l) = 0.5 * (elasticity_tensor(i, j, k, l) - elasticity_tensor(j, i, k, l));
     141             : 
     142         128 :   cto += anti;
     143             : }