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 "ComputeMultipleInelasticStress.h"
      11             : 
      12             : #include "StressUpdateBase.h"
      13             : #include "MooseException.h"
      14             : #include "DamageBase.h"
      15             : #include "libmesh/int_range.h"
      16             : 
      17             : registerMooseObject("SolidMechanicsApp", ComputeMultipleInelasticStress);
      18             : 
      19             : InputParameters
      20        5858 : ComputeMultipleInelasticStress::validParams()
      21             : {
      22        5858 :   InputParameters params = ComputeMultipleInelasticStressBase::validParams();
      23       11716 :   params.addRequiredParam<std::vector<MaterialName>>(
      24             :       "inelastic_models",
      25             :       "The material objects to use to calculate stress and inelastic strains. "
      26             :       "Note: specify creep models first and plasticity models second.");
      27        5858 :   return params;
      28           0 : }
      29             : 
      30        4388 : ComputeMultipleInelasticStress::ComputeMultipleInelasticStress(const InputParameters & parameters)
      31        4388 :   : ComputeMultipleInelasticStressBase(parameters)
      32             : {
      33        4388 : }
      34             : 
      35             : std::vector<MaterialName>
      36        4194 : ComputeMultipleInelasticStress::getInelasticModelNames()
      37             : {
      38       12582 :   return getParam<std::vector<MaterialName>>("inelastic_models");
      39             : }
      40             : 
      41             : void
      42     1526416 : ComputeMultipleInelasticStress::updateQpState(RankTwoTensor & elastic_strain_increment,
      43             :                                               RankTwoTensor & combined_inelastic_strain_increment)
      44             : {
      45     1526416 :   if (_internal_solve_full_iteration_history == true)
      46             :   {
      47           0 :     _console << std::endl
      48           0 :              << "iteration output for ComputeMultipleInelasticStress solve:"
      49           0 :              << " time=" << _t << " int_pt=" << _qp << std::endl;
      50             :   }
      51             :   Real l2norm_delta_stress;
      52             :   Real first_l2norm_delta_stress = 1.0;
      53     1526416 :   unsigned int counter = 0;
      54             : 
      55             :   std::vector<RankTwoTensor> inelastic_strain_increment;
      56     1526416 :   inelastic_strain_increment.resize(_num_models);
      57             : 
      58     4579248 :   for (const auto i_rmm : index_range(_models))
      59             :     inelastic_strain_increment[i_rmm].zero();
      60             : 
      61     1526416 :   RankTwoTensor stress_max, stress_min;
      62             : 
      63             :   do
      64             :   {
      65    18584148 :     for (const auto i_rmm : index_range(_models))
      66             :     {
      67    12389432 :       _models[i_rmm]->setQp(_qp);
      68             : 
      69             :       // initially assume the strain is completely elastic
      70    12389432 :       elastic_strain_increment = _strain_increment[_qp];
      71             :       // and subtract off all inelastic strain increments calculated so far
      72             :       // except the one that we're about to calculate
      73    37168296 :       for (const auto j_rmm : make_range(_num_models))
      74    24778864 :         if (i_rmm != j_rmm)
      75    12389432 :           elastic_strain_increment -= inelastic_strain_increment[j_rmm];
      76             : 
      77             :       // form the trial stress, with the check for changed elasticity constants
      78    12389432 :       if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
      79    12383672 :         _stress[_qp] =
      80    12383672 :             _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
      81             :       else
      82             :       {
      83        5760 :         if (_damage_model)
      84           0 :           _stress[_qp] = _undamaged_stress_old + _elasticity_tensor[_qp] * elastic_strain_increment;
      85             :         else
      86        5760 :           _stress[_qp] = _stress_old[_qp] + _elasticity_tensor[_qp] * elastic_strain_increment;
      87             :       }
      88             : 
      89             :       // given a trial stress (_stress[_qp]) and a strain increment (elastic_strain_increment)
      90             :       // let the i^th model produce an admissible stress (as _stress[_qp]), and decompose
      91             :       // the strain increment into an elastic part (elastic_strain_increment) and an
      92             :       // inelastic part (inelastic_strain_increment[i_rmm])
      93    12389432 :       computeAdmissibleState(i_rmm,
      94             :                              elastic_strain_increment,
      95             :                              inelastic_strain_increment[i_rmm],
      96             :                              _consistent_tangent_operator[i_rmm]);
      97             : 
      98    12389432 :       if (i_rmm == 0)
      99             :       {
     100     6194716 :         stress_max = _stress[_qp];
     101     6194716 :         stress_min = _stress[_qp];
     102             :       }
     103             :       else
     104             :       {
     105    24778864 :         for (const auto i : make_range(Moose::dim))
     106    74336592 :           for (const auto j : make_range(Moose::dim))
     107    55752444 :             if (_stress[_qp](i, j) > stress_max(i, j))
     108    23755300 :               stress_max(i, j) = _stress[_qp](i, j);
     109    31997144 :             else if (stress_min(i, j) > _stress[_qp](i, j))
     110    18819504 :               stress_min(i, j) = _stress[_qp](i, j);
     111             :       }
     112             :     }
     113             : 
     114             :     // now check convergence in the stress:
     115             :     // once the change in stress is within tolerance after each recompute material
     116             :     // consider the stress to be converged
     117     6194716 :     l2norm_delta_stress = (stress_max - stress_min).L2norm();
     118     6194716 :     if (counter == 0 && l2norm_delta_stress > 0.0)
     119             :       first_l2norm_delta_stress = l2norm_delta_stress;
     120             : 
     121     6194716 :     if (_internal_solve_full_iteration_history == true)
     122             :     {
     123           0 :       _console << "stress iteration number = " << counter << "\n"
     124           0 :                << " relative l2 norm delta stress = "
     125           0 :                << (0 == first_l2norm_delta_stress ? 0
     126           0 :                                                   : l2norm_delta_stress / first_l2norm_delta_stress)
     127           0 :                << "\n"
     128           0 :                << " stress convergence relative tolerance = " << _relative_tolerance << "\n"
     129           0 :                << " absolute l2 norm delta stress = " << l2norm_delta_stress << "\n"
     130           0 :                << " stress convergence absolute tolerance = " << _absolute_tolerance << std::endl;
     131             :     }
     132     6194716 :     ++counter;
     133     6191836 :   } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
     134    10877400 :            (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
     135     4668300 :            _num_models != 1);
     136             : 
     137     1526416 :   if (counter == _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
     138           0 :       (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance)
     139           0 :     throw MooseException("Max stress iteration hit during ComputeMultipleInelasticStress solve!");
     140             : 
     141             :   combined_inelastic_strain_increment.zero();
     142     4579248 :   for (const auto i_rmm : make_range(_num_models))
     143             :     combined_inelastic_strain_increment +=
     144     3052832 :         _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
     145             : 
     146     1526416 :   if (_fe_problem.currentlyComputingJacobian())
     147      139104 :     computeQpJacobianMult();
     148             : 
     149     1526416 :   _material_timestep_limit[_qp] = 0.0;
     150     4579248 :   for (const auto i_rmm : make_range(_num_models))
     151     3052832 :     _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
     152             : 
     153     1526416 :   if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
     154      116594 :     _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
     155             :   else
     156     1409822 :     _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
     157     1526416 : }