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 "ComputeCreepPlasticityStress.h"
      11             : 
      12             : #include "ElasticityTensorTools.h"
      13             : 
      14             : #include "MooseException.h"
      15             : #include "libmesh/int_range.h"
      16             : 
      17             : registerMooseObject("SolidMechanicsApp", ComputeCreepPlasticityStress);
      18             : 
      19             : InputParameters
      20          48 : ComputeCreepPlasticityStress::validParams()
      21             : {
      22          48 :   InputParameters params = ComputeMultipleInelasticStressBase::validParams();
      23          48 :   params.addClassDescription("Compute state (stress and internal parameters such as inelastic "
      24             :                              "strains and internal parameters) using an Newton process for one "
      25             :                              "creep and one plasticity model");
      26          96 :   params.addRequiredParam<MaterialName>("creep_model",
      27             :                                         "Creep model that derives from PowerLawCreepStressUpdate.");
      28          96 :   params.addRequiredParam<MaterialName>(
      29             :       "plasticity_model", "Plasticity model that derives from IsotropicPlasticityStressUpdate.");
      30             : 
      31          48 :   return params;
      32           0 : }
      33             : 
      34          36 : ComputeCreepPlasticityStress::ComputeCreepPlasticityStress(const InputParameters & parameters)
      35          36 :   : ComputeMultipleInelasticStressBase(parameters)
      36             : {
      37          36 : }
      38             : 
      39             : std::vector<MaterialName>
      40          36 : ComputeCreepPlasticityStress::getInelasticModelNames()
      41             : {
      42             :   std::vector<MaterialName> names = {getParam<MaterialName>("creep_model"),
      43         108 :                                      getParam<MaterialName>("plasticity_model")};
      44          36 :   return names;
      45         108 : }
      46             : 
      47             : void
      48          36 : ComputeCreepPlasticityStress::initialSetup()
      49             : {
      50          36 :   ComputeMultipleInelasticStressBase::initialSetup();
      51             : 
      52          36 :   if (_models.size() != 2)
      53           0 :     mooseError("Error in ComputeCreepPlasticityStress: two models are required.");
      54             : 
      55          36 :   _creep_model = dynamic_cast<PowerLawCreepStressUpdate *>(_models[0]);
      56          36 :   if (!_creep_model)
      57           0 :     mooseError("Model " + getParam<MaterialName>("creep_model") +
      58             :                " is not a compatible creep model in ComputeCreepPlasticityStress.");
      59             : 
      60          36 :   _plasticity_model = dynamic_cast<IsotropicPlasticityStressUpdate *>(_models[1]);
      61          36 :   if (!_plasticity_model)
      62           0 :     mooseError("Model " + getParam<MaterialName>("plasticity_model") +
      63             :                " is not a compatible plasticity model in ComputeCreepPlasticityStress.");
      64          36 : }
      65             : 
      66             : void
      67      126240 : ComputeCreepPlasticityStress::updateQpState(RankTwoTensor & elastic_strain_increment,
      68             :                                             RankTwoTensor & combined_inelastic_strain_increment)
      69             : {
      70      126240 :   if (_internal_solve_full_iteration_history == true)
      71             :   {
      72           0 :     _console << std::endl
      73           0 :              << "iteration output for ComputeCreepPlasticityStress solve:"
      74           0 :              << " time=" << _t << " int_pt=" << _qp << std::endl;
      75             :   }
      76             : 
      77      378720 :   for (const auto i_rmm : index_range(_models))
      78             :     _inelastic_strain_increment[i_rmm].zero();
      79             : 
      80      126240 :   elastic_strain_increment = _strain_increment[_qp];
      81      126240 :   computeStress(elastic_strain_increment, _stress[_qp]);
      82             : 
      83      126240 :   const RankTwoTensor trial_stress = _stress[_qp];
      84      126240 :   const RankTwoTensor deviatoric_trial_stress = trial_stress.deviatoric();
      85             :   const Real effective_trial_stress =
      86      126240 :       std::sqrt(1.5 * deviatoric_trial_stress.doubleContraction(deviatoric_trial_stress));
      87             : 
      88             :   //
      89             :   // First, compute creep response assuming no plasticity
      90             :   //
      91      126240 :   _creep_model->setQp(_qp);
      92      126240 :   computeAdmissibleState(
      93             :       0, elastic_strain_increment, _inelastic_strain_increment[0], _consistent_tangent_operator[0]);
      94      126240 :   _effective_creep_strain_increment = _creep_model->effectiveInelasticStrainIncrement();
      95             : 
      96             :   //
      97             :   // Now check the plasticity model.
      98             :   // If no yielding, we are done.
      99             :   //
     100      126240 :   _plasticity_model->setQp(_qp);
     101      126240 :   computeAdmissibleState(
     102             :       1, elastic_strain_increment, _inelastic_strain_increment[1], _consistent_tangent_operator[1]);
     103      126240 :   _effective_plastic_strain_increment = _plasticity_model->effectiveInelasticStrainIncrement();
     104             :   const Real yield_condition = _plasticity_model->yieldCondition();
     105             : 
     106      126240 :   if (yield_condition > 0.0) // yielding even with maximum creep strain
     107             :   {
     108       45608 :     const Real max_effective_creep_strain_increment = _effective_creep_strain_increment;
     109             : 
     110             :     //
     111             :     // Compute plastic response given no creep
     112             :     //
     113       45608 :     elastic_strain_increment = _strain_increment[_qp];
     114       45608 :     _stress[_qp] = trial_stress;
     115       45608 :     computeAdmissibleState(1,
     116             :                            elastic_strain_increment,
     117             :                            _inelastic_strain_increment[1],
     118             :                            _consistent_tangent_operator[1]);
     119             :     const Real max_effective_plastic_strain_increment =
     120       45608 :         _plasticity_model->effectiveInelasticStrainIncrement();
     121             : 
     122             :     //
     123             :     // Reset plasticity model
     124             :     //
     125       45608 :     _plasticity_model->computeStressInitialize(effective_trial_stress, _elasticity_tensor[_qp]);
     126             : 
     127             :     //
     128             :     // The residual for the creep law is (creep rate)*dt - (creep strain increment)
     129             :     // We want the creep rate calculation to depend on both creep and plasticity.
     130             :     // Since we send in the total scalar inelastic strain (creep and plasticity), we need to correct
     131             :     // the second term by adding the plastic strain.
     132             :     //
     133       45608 :     Real creep_residual = _creep_model->computeResidual(effective_trial_stress,
     134       91216 :                                                         _effective_creep_strain_increment +
     135       45608 :                                                             _effective_plastic_strain_increment);
     136       45608 :     creep_residual += _effective_plastic_strain_increment;
     137             : 
     138             :     //
     139             :     // We want the plasticity resdiual to be (effective_trial_stress - r - yield_stress)/3G - (total
     140             :     // inelastic increment)
     141             :     // The standard residual for plasticity is (effective_trial_stress - r - yield_stress)/3G -
     142             :     // (plasticity strain increment)
     143             :     // Since we send in the plastic inelastic strain (for calculation of r), we must subtract the
     144             :     // creep strain to correct the final term in our desired residual
     145             :     //
     146       45608 :     Real plasticity_residual = _plasticity_model->computeResidual(
     147       45608 :         effective_trial_stress, _effective_plastic_strain_increment);
     148       45608 :     plasticity_residual -= _effective_creep_strain_increment;
     149             : 
     150       45608 :     unsigned int counter = 0;
     151       45608 :     Real residual = 0;
     152             :     const Real threeG =
     153       45608 :         3.0 * ElasticityTensorTools::getIsotropicShearModulus(_elasticity_tensor[_qp]);
     154       45608 :     const Real reference_residual = effective_trial_stress / threeG;
     155             : 
     156             :     do
     157             :     {
     158             :       //
     159             :       // Solve Ax=b where x is the vector of creep and plastic inelastic strains
     160             :       //
     161             :       // x1 => creep
     162             :       // x2 => plasticity
     163             :       //
     164             :       //
     165             :       // A11 (creep,creep)
     166             :       //
     167       77560 :       Real A11 = _creep_model->computeDerivative(effective_trial_stress,
     168      155120 :                                                  _effective_creep_strain_increment +
     169       77560 :                                                      _effective_plastic_strain_increment);
     170             :       //
     171             :       // A12 (creep,plasticity)
     172             :       //
     173       77560 :       Real A12 = A11 + 1.0;
     174             :       //
     175             :       // A21 (plasticity,creep)
     176             :       //
     177             :       Real A21 = -1.0;
     178             :       //
     179             :       // A22 (plasticity,plasticity)
     180             :       //
     181       77560 :       Real A22 = _plasticity_model->computeDerivative(effective_trial_stress,
     182             :                                                       _effective_plastic_strain_increment);
     183             : 
     184       77560 :       Real rhs1 = creep_residual;
     185       77560 :       Real rhs2 = plasticity_residual;
     186             : 
     187             :       //
     188             :       // Solve
     189             :       //
     190             :       // A11*a + A21 = 0
     191             :       // A11*a = -A21
     192             :       // a = -A21/A11
     193       77560 :       const Real a = -A21 / A11;
     194             :       A21 = 0;
     195       77560 :       A22 += a * A12;
     196       77560 :       rhs2 += a * rhs1;
     197             : 
     198       77560 :       Real x2 = rhs2 / A22;
     199       77560 :       Real x1 = (rhs1 - A12 * x2) / A11;
     200             : 
     201       77560 :       while (_effective_creep_strain_increment - x1 < 0)
     202           0 :         x1 *= 0.5;
     203       77560 :       while (_effective_plastic_strain_increment - x2 < 0)
     204           0 :         x2 *= 0.5;
     205             : 
     206       77560 :       while (_effective_creep_strain_increment - x1 > max_effective_creep_strain_increment)
     207           0 :         x1 *= 0.5;
     208       77560 :       while (_effective_plastic_strain_increment - x2 > max_effective_plastic_strain_increment)
     209           0 :         x2 *= 0.5;
     210             : 
     211             :       // This check is to avoid a fpe in the creep law.
     212             :       // Maybe it is better to check for this condition in the creep law
     213       77560 :       if (effective_trial_stress < threeG * (_effective_creep_strain_increment - x1 +
     214       77560 :                                              _effective_plastic_strain_increment - x2))
     215             :       {
     216             :         const Real sum =
     217             :             _effective_creep_strain_increment - x1 + _effective_plastic_strain_increment - x2;
     218           0 :         const Real factor = 0.9 * effective_trial_stress / (threeG * sum);
     219           0 :         const Real cNew = (_effective_creep_strain_increment - x1) * factor;
     220           0 :         const Real pNew = (_effective_plastic_strain_increment - x2) * factor;
     221           0 :         x1 = _effective_creep_strain_increment - cNew;
     222           0 :         x2 = _effective_plastic_strain_increment - pNew;
     223             :       }
     224             : 
     225       77560 :       _effective_creep_strain_increment -= x1;
     226       77560 :       _effective_plastic_strain_increment -= x2;
     227             : 
     228             :       //
     229             :       // The residual for the creep law is (creep rate)*dt - (creep strain increment)
     230             :       //
     231       77560 :       creep_residual = _creep_model->computeResidual(effective_trial_stress,
     232       77560 :                                                      _effective_creep_strain_increment +
     233             :                                                          _effective_plastic_strain_increment);
     234       77560 :       creep_residual += _effective_plastic_strain_increment;
     235             :       //
     236             :       // The residual for plasticity is (effective_trial_stress - r - yield_stress)/3G - scalar
     237             :       //
     238       77560 :       plasticity_residual = _plasticity_model->computeResidual(effective_trial_stress,
     239             :                                                                _effective_plastic_strain_increment);
     240       77560 :       plasticity_residual -= _effective_creep_strain_increment;
     241             : 
     242       77560 :       residual =
     243       77560 :           std::sqrt(creep_residual * creep_residual + plasticity_residual * plasticity_residual);
     244             : 
     245       77560 :       if (_internal_solve_full_iteration_history == true)
     246             :       {
     247           0 :         _console << "stress iteration number = " << counter << "\n relative residual = "
     248           0 :                  << (0 == reference_residual ? 0 : residual / reference_residual)
     249           0 :                  << "\n stress convergence relative tolerance = " << _relative_tolerance
     250           0 :                  << "\n absolute residual = " << residual
     251           0 :                  << "\n creep residual = " << creep_residual
     252           0 :                  << "\n plasticity residual = " << plasticity_residual
     253           0 :                  << "\n creep iteration increment = " << x1
     254           0 :                  << "\n plasticity iteration increment = " << x2
     255           0 :                  << "\n creep scalar strain = " << _effective_creep_strain_increment
     256           0 :                  << "\n plasticity scalar strain = " << _effective_plastic_strain_increment
     257           0 :                  << "\n max creep scalar strain = " << max_effective_creep_strain_increment
     258           0 :                  << "\n max plasticity scalar strain = " << max_effective_plastic_strain_increment
     259           0 :                  << "\n stress convergence absolute tolerance = " << _absolute_tolerance
     260           0 :                  << std::endl;
     261             :       }
     262       77560 :       ++counter;
     263       77560 :     } while (counter < _max_iterations && residual > _absolute_tolerance &&
     264       31952 :              (residual / reference_residual) > _relative_tolerance);
     265             : 
     266       45608 :     if (counter == _max_iterations && residual > _absolute_tolerance &&
     267           0 :         (residual / reference_residual) > _relative_tolerance)
     268             :     {
     269           0 :       std::stringstream msg;
     270           0 :       msg << "\n relative residual = "
     271           0 :           << (0 == reference_residual ? 0 : residual / reference_residual)
     272           0 :           << "\n stress convergence relative tolerance = " << _relative_tolerance
     273           0 :           << "\n absolute residual = " << residual << "\n creep residual = " << creep_residual
     274           0 :           << "\n plasticity residual = " << plasticity_residual
     275           0 :           << "\n creep scalar strain = " << _effective_creep_strain_increment
     276           0 :           << "\n plasticity scalar strain = " << _effective_plastic_strain_increment
     277           0 :           << "\n max creep scalar strain = " << max_effective_creep_strain_increment
     278           0 :           << "\n max plasticity scalar strain = " << max_effective_plastic_strain_increment
     279           0 :           << "\n stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
     280             : 
     281           0 :       throw MooseException("Max stress iteration hit during ComputeCreepPlasticityStress solve! " +
     282           0 :                            _name + msg.str());
     283           0 :     }
     284             :   }
     285      126240 :   computeInelasticStrainIncrements(effective_trial_stress, deviatoric_trial_stress);
     286             : 
     287      126240 :   finalizeConstitutiveModels();
     288             : 
     289             :   combined_inelastic_strain_increment.zero();
     290      378720 :   for (const auto i_rmm : make_range(_num_models))
     291      252480 :     combined_inelastic_strain_increment += _inelastic_strain_increment[i_rmm];
     292             : 
     293      126240 :   if (yield_condition > 0.0)
     294             :   {
     295       45608 :     elastic_strain_increment = _strain_increment[_qp] - combined_inelastic_strain_increment;
     296       45608 :     _stress[_qp] = _elasticity_tensor[_qp] * (elastic_strain_increment + _elastic_strain_old[_qp]);
     297             : 
     298       45608 :     computeTangentOperators();
     299             :   }
     300             :   else
     301             :   {
     302             :     // The tangent for creep was called in the initial creep solve.
     303             : 
     304             :     // Set the tangent for plasticity to zero
     305       80632 :     _consistent_tangent_operator[1].zero();
     306             :   }
     307             : 
     308      126240 :   if (_fe_problem.currentlyComputingJacobian())
     309       18896 :     computeQpJacobianMult();
     310             : 
     311      126240 :   _material_timestep_limit[_qp] = 0.0;
     312      378720 :   for (const auto i_rmm : make_range(_num_models))
     313      252480 :     _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
     314             : 
     315      126240 :   if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
     316        1632 :     _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
     317             :   else
     318      124608 :     _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
     319      126240 : }
     320             : 
     321             : void
     322      217456 : ComputeCreepPlasticityStress::computeStress(const RankTwoTensor & elastic_strain_increment,
     323             :                                             RankTwoTensor & stress)
     324             : {
     325             :   // form the stress, with the check for changed elasticity constants
     326      217456 :   if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
     327      217456 :     stress = _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
     328             :   else
     329             :   {
     330           0 :     if (_damage_model)
     331           0 :       stress = _undamaged_stress_old + _elasticity_tensor[_qp] * elastic_strain_increment;
     332             :     else
     333           0 :       stress = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
     334             :   }
     335      217456 : }
     336             : 
     337             : void
     338       45608 : ComputeCreepPlasticityStress::computeTangentOperators()
     339             : {
     340       45608 :   auto elastic_strain_inc = _strain_increment[_qp] - _inelastic_strain_increment[1];
     341             : 
     342       45608 :   RankTwoTensor stress;
     343       45608 :   computeStress(elastic_strain_inc, stress);
     344       45608 :   RankTwoTensor deviatoric_stress = stress.deviatoric();
     345       45608 :   Real effective_stress = std::sqrt(1.5 * deviatoric_stress.doubleContraction(deviatoric_stress));
     346             : 
     347       45608 :   _creep_model->computeTangentOperator(
     348       45608 :       effective_stress, _stress[_qp], _consistent_tangent_operator[0]);
     349             : 
     350       45608 :   elastic_strain_inc = _strain_increment[_qp] - _inelastic_strain_increment[0];
     351             : 
     352       45608 :   computeStress(elastic_strain_inc, stress);
     353       45608 :   deviatoric_stress = stress.deviatoric();
     354       45608 :   effective_stress = std::sqrt(1.5 * deviatoric_stress.doubleContraction(deviatoric_stress));
     355             : 
     356       45608 :   _plasticity_model->computeTangentOperator(
     357       45608 :       effective_stress, _stress[_qp], _consistent_tangent_operator[1]);
     358       45608 : }
     359             : 
     360             : void
     361      126240 : ComputeCreepPlasticityStress::computeInelasticStrainIncrements(
     362             :     Real effective_trial_stress, const RankTwoTensor & deviatoric_trial_stress)
     363             : {
     364      126240 :   if (!MooseUtils::absoluteFuzzyEqual(effective_trial_stress, 0.0))
     365             :   {
     366      126144 :     if (_effective_creep_strain_increment != 0.0)
     367      111624 :       _inelastic_strain_increment[0] =
     368      111624 :           deviatoric_trial_stress *
     369      111624 :           (1.5 * _effective_creep_strain_increment / effective_trial_stress);
     370             :     else
     371             :       _inelastic_strain_increment[0].zero();
     372             : 
     373      126144 :     if (_effective_plastic_strain_increment != 0.0)
     374       45608 :       _inelastic_strain_increment[1] =
     375       45608 :           deviatoric_trial_stress *
     376       45608 :           (1.5 * _effective_plastic_strain_increment / effective_trial_stress);
     377             :     else
     378             :       _inelastic_strain_increment[1].zero();
     379             :   }
     380             :   else
     381             :   {
     382             :     _inelastic_strain_increment[0].zero();
     383             :     _inelastic_strain_increment[1].zero();
     384             :   }
     385      126240 : }
     386             : 
     387             : void
     388      126240 : ComputeCreepPlasticityStress::finalizeConstitutiveModels()
     389             : {
     390      126240 :   _creep_model->updateEffectiveInelasticStrainIncrement(_effective_creep_strain_increment);
     391      126240 :   _plasticity_model->updateEffectiveInelasticStrainIncrement(_effective_plastic_strain_increment);
     392      126240 :   _creep_model->updateEffectiveInelasticStrain(_effective_creep_strain_increment);
     393      126240 :   _plasticity_model->updateEffectiveInelasticStrain(_effective_plastic_strain_increment);
     394      126240 :   _creep_model->resetIncrementalMaterialProperties();
     395      126240 :   _creep_model->computeStressFinalize(_inelastic_strain_increment[0]);
     396      126240 :   _plasticity_model->computeStressFinalize(_inelastic_strain_increment[1]);
     397      126240 : }