Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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("TensorMechanicsApp", ComputeMultipleInelasticStress);
      18             : 
      19             : InputParameters
      20        2761 : ComputeMultipleInelasticStress::validParams()
      21             : {
      22        2761 :   InputParameters params = ComputeMultipleInelasticStressBase::validParams();
      23        5522 :   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        2761 :   return params;
      28           0 : }
      29             : 
      30        2068 : ComputeMultipleInelasticStress::ComputeMultipleInelasticStress(const InputParameters & parameters)
      31        2068 :   : ComputeMultipleInelasticStressBase(parameters)
      32             : {
      33        2068 : }
      34             : 
      35             : std::vector<MaterialName>
      36        1980 : ComputeMultipleInelasticStress::getInelasticModelNames()
      37             : {
      38        5940 :   return getParam<std::vector<MaterialName>>("inelastic_models");
      39             : }
      40             : 
      41             : void
      42      294168 : ComputeMultipleInelasticStress::updateQpState(RankTwoTensor & elastic_strain_increment,
      43             :                                               RankTwoTensor & combined_inelastic_strain_increment)
      44             : {
      45      294168 :   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      294168 :   unsigned int counter = 0;
      54             : 
      55             :   std::vector<RankTwoTensor> inelastic_strain_increment;
      56      294168 :   inelastic_strain_increment.resize(_num_models);
      57             : 
      58      882504 :   for (const auto i_rmm : index_range(_models))
      59             :     inelastic_strain_increment[i_rmm].zero();
      60             : 
      61      294168 :   RankTwoTensor stress_max, stress_min;
      62             : 
      63             :   do
      64             :   {
      65     7924842 :     for (const auto i_rmm : index_range(_models))
      66             :     {
      67     5283228 :       _models[i_rmm]->setQp(_qp);
      68             : 
      69             :       // initially assume the strain is completely elastic
      70     5283228 :       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    15849684 :       for (const auto j_rmm : make_range(_num_models))
      74    10566456 :         if (i_rmm != j_rmm)
      75     5283228 :           elastic_strain_increment -= inelastic_strain_increment[j_rmm];
      76             : 
      77             :       // form the trial stress, with the check for changed elasticity constants
      78     5283228 :       if (_is_elasticity_tensor_guaranteed_isotropic || !_perform_finite_strain_rotations)
      79     5280284 :         _stress[_qp] =
      80     5280284 :             _elasticity_tensor[_qp] * (_elastic_strain_old[_qp] + elastic_strain_increment);
      81             :       else
      82             :       {
      83        2944 :         if (_damage_model)
      84           0 :           _stress[_qp] = _undamaged_stress_old + _elasticity_tensor[_qp] * elastic_strain_increment;
      85             :         else
      86        2944 :           _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     5283228 :       computeAdmissibleState(i_rmm,
      94             :                              elastic_strain_increment,
      95             :                              inelastic_strain_increment[i_rmm],
      96             :                              _consistent_tangent_operator[i_rmm]);
      97             : 
      98     5283228 :       if (i_rmm == 0)
      99             :       {
     100     2641614 :         stress_max = _stress[_qp];
     101     2641614 :         stress_min = _stress[_qp];
     102             :       }
     103             :       else
     104             :       {
     105    10566456 :         for (const auto i : make_range(Moose::dim))
     106    31699368 :           for (const auto j : make_range(Moose::dim))
     107    23774526 :             if (_stress[_qp](i, j) > stress_max(i, j))
     108    11795874 :               stress_max(i, j) = _stress[_qp](i, j);
     109    11978652 :             else if (stress_min(i, j) > _stress[_qp](i, j))
     110     9605136 :               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     2641614 :     l2norm_delta_stress = (stress_max - stress_min).L2norm();
     118     2641614 :     if (counter == 0 && l2norm_delta_stress > 0.0)
     119             :       first_l2norm_delta_stress = l2norm_delta_stress;
     120             : 
     121     2641614 :     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     2641614 :     ++counter;
     133     2640142 :   } while (counter < _max_iterations && l2norm_delta_stress > _absolute_tolerance &&
     134     4996588 :            (l2norm_delta_stress / first_l2norm_delta_stress) > _relative_tolerance &&
     135     2347446 :            _num_models != 1);
     136             : 
     137      294168 :   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      882504 :   for (const auto i_rmm : make_range(_num_models))
     143             :     combined_inelastic_strain_increment +=
     144      588336 :         _inelastic_weights[i_rmm] * inelastic_strain_increment[i_rmm];
     145             : 
     146      294168 :   if (_fe_problem.currentlyComputingJacobian())
     147       47448 :     computeQpJacobianMult();
     148             : 
     149      294168 :   _material_timestep_limit[_qp] = 0.0;
     150      882504 :   for (const auto i_rmm : make_range(_num_models))
     151      588336 :     _material_timestep_limit[_qp] += 1.0 / _models[i_rmm]->computeTimeStepLimit();
     152             : 
     153      294168 :   if (MooseUtils::absoluteFuzzyEqual(_material_timestep_limit[_qp], 0.0))
     154       64320 :     _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
     155             :   else
     156      229848 :     _material_timestep_limit[_qp] = 1.0 / _material_timestep_limit[_qp];
     157      294168 : }