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 "ComputeIsotropicElasticityTensor.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", ComputeIsotropicElasticityTensor);
      13             : registerMooseObject("SolidMechanicsApp", ADComputeIsotropicElasticityTensor);
      14             : registerMooseObject("SolidMechanicsApp", SymmetricIsotropicElasticityTensor);
      15             : registerMooseObject("SolidMechanicsApp", ADSymmetricIsotropicElasticityTensor);
      16             : 
      17             : template <bool is_ad, typename T>
      18             : InputParameters
      19       22422 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::validParams()
      20             : {
      21       22422 :   InputParameters params = ComputeElasticityTensorBase::validParams();
      22       22422 :   params.addClassDescription("Compute a constant isotropic elasticity tensor.");
      23       44844 :   params.addParam<Real>("bulk_modulus", -1, "The bulk modulus for the material.");
      24       44844 :   params.addParam<Real>("lambda", -1, "Lame's first constant for the material.");
      25       44844 :   params.addParam<Real>("poissons_ratio", -1, "Poisson's ratio for the material.");
      26       44844 :   params.addParam<Real>("shear_modulus", -1, "The shear modulus of the material.");
      27       44844 :   params.addParam<Real>("youngs_modulus", -1, "Young's modulus of the material.");
      28       44844 :   params.declareControllable("bulk_modulus lambda poissons_ratio shear_modulus youngs_modulus");
      29       22422 :   return params;
      30           0 : }
      31             : 
      32             : template <bool is_ad, typename T>
      33       16804 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::ComputeIsotropicElasticityTensorTempl(
      34             :     const InputParameters & parameters)
      35             :   : ComputeElasticityTensorBaseTempl<is_ad, T>(parameters),
      36       16804 :     _bulk_modulus_set(parameters.isParamSetByUser("bulk_modulus")),
      37       16804 :     _lambda_set(parameters.isParamSetByUser("lambda")),
      38       16804 :     _poissons_ratio_set(parameters.isParamSetByUser("poissons_ratio")),
      39       16804 :     _shear_modulus_set(parameters.isParamSetByUser("shear_modulus")),
      40       16804 :     _youngs_modulus_set(parameters.isParamSetByUser("youngs_modulus")),
      41       33608 :     _bulk_modulus(this->template getParam<Real>("bulk_modulus")),
      42       33608 :     _lambda(this->template getParam<Real>("lambda")),
      43       33608 :     _poissons_ratio(this->template getParam<Real>("poissons_ratio")),
      44       33608 :     _shear_modulus(this->template getParam<Real>("shear_modulus")),
      45       33608 :     _youngs_modulus(this->template getParam<Real>("youngs_modulus")),
      46       33608 :     _effective_stiffness_local(parameters.isParamValid("effective_stiffness_local"))
      47             : {
      48       16804 :   unsigned int num_elastic_constants = _bulk_modulus_set + _lambda_set + _poissons_ratio_set +
      49       16804 :                                        _shear_modulus_set + _youngs_modulus_set;
      50       16804 :   if (num_elastic_constants != 2)
      51           4 :     mooseError("Exactly two elastic constants must be defined for material '" + name() + "'.");
      52             : 
      53             :   // all tensors created by this class are always isotropic
      54       33604 :   issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);
      55       16802 :   issueGuarantee("effective_stiffness", Guarantee::ISOTROPIC);
      56       33604 :   if (!isParamValid("elasticity_tensor_prefactor"))
      57       33388 :     issueGuarantee(_elasticity_tensor_name, Guarantee::CONSTANT_IN_TIME);
      58             : 
      59       16802 :   if (_bulk_modulus_set && _bulk_modulus <= 0.0)
      60           4 :     mooseError("Bulk modulus must be positive in material '" + name() + "'.");
      61             : 
      62       16800 :   if (_poissons_ratio_set && (_poissons_ratio <= -1.0 || _poissons_ratio >= 0.5))
      63           4 :     mooseError("Poissons ratio must be greater than -1 and less than 0.5 in "
      64             :                "material '" +
      65             :                name() + "'.");
      66             : 
      67       16798 :   if (_shear_modulus_set && _shear_modulus < 0.0)
      68           4 :     mooseError("Shear modulus must not be negative in material '" + name() + "'.");
      69             : 
      70       16796 :   if (_youngs_modulus_set && _youngs_modulus <= 0.0)
      71           4 :     mooseError("Youngs modulus must be positive in material '" + name() + "'.");
      72       16794 : }
      73             : 
      74             : template <bool is_ad, typename T>
      75             : void
      76     3586644 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::residualSetup()
      77             : {
      78     3586644 :   std::vector<Real> iso_const(2);
      79             :   Real elas_mod;
      80             :   Real poiss_rat;
      81             : 
      82     3586644 :   if (_youngs_modulus_set && _poissons_ratio_set)
      83             :   {
      84     3439626 :     _Cijkl.fillSymmetricIsotropicEandNu(_youngs_modulus, _poissons_ratio);
      85     3439626 :     _effective_stiffness_local =
      86     6879252 :         std::max(std::sqrt((_youngs_modulus * (1 - _poissons_ratio)) /
      87     3439626 :                            ((1 + _poissons_ratio) * (1 - 2 * _poissons_ratio))),
      88     3439626 :                  std::sqrt(_youngs_modulus / (2 * (1 + _poissons_ratio))));
      89             :     return;
      90             :   }
      91             : 
      92      147018 :   if (_lambda_set && _shear_modulus_set)
      93             :   {
      94      112212 :     iso_const[0] = _lambda;
      95      112212 :     iso_const[1] = _shear_modulus;
      96      112212 :     elas_mod = (_shear_modulus * (3 * _lambda + 2 * _shear_modulus)) / (_lambda + _shear_modulus);
      97      112212 :     poiss_rat = _lambda / (2 * (_lambda + _shear_modulus));
      98      112212 :     _effective_stiffness_local =
      99      224424 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     100      112212 :                  std::sqrt(_shear_modulus));
     101             :   }
     102       34806 :   else if (_shear_modulus_set && _bulk_modulus_set)
     103             :   {
     104         630 :     iso_const[0] = _bulk_modulus - 2.0 / 3.0 * _shear_modulus;
     105         630 :     iso_const[1] = _shear_modulus;
     106         630 :     elas_mod = (9 * _bulk_modulus * _shear_modulus) / (3 * _bulk_modulus + _shear_modulus);
     107         630 :     poiss_rat =
     108         630 :         (3 * _bulk_modulus - 2 * _shear_modulus) / (2 * (3 * _bulk_modulus + _shear_modulus));
     109         630 :     _effective_stiffness_local =
     110        1260 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     111         630 :                  std::sqrt(_shear_modulus));
     112             :   }
     113       34176 :   else if (_poissons_ratio_set && _bulk_modulus_set)
     114             :   {
     115           0 :     iso_const[0] = 3.0 * _bulk_modulus * _poissons_ratio / (1.0 + _poissons_ratio);
     116           0 :     iso_const[1] =
     117           0 :         3.0 * _bulk_modulus * (1.0 - 2.0 * _poissons_ratio) / (2.0 * (1.0 + _poissons_ratio));
     118           0 :     elas_mod = 3 * _bulk_modulus * (1 - 2 * _poissons_ratio);
     119             :     poiss_rat = _poissons_ratio;
     120           0 :     _effective_stiffness_local =
     121           0 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     122           0 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     123             :   }
     124       34176 :   else if (_lambda_set && _bulk_modulus_set)
     125             :   {
     126           0 :     iso_const[0] = _lambda;
     127           0 :     iso_const[1] = 3.0 * (_bulk_modulus - _lambda) / 2.0;
     128           0 :     elas_mod = (9 * _bulk_modulus * (_bulk_modulus - _lambda)) / (3 * _bulk_modulus - _lambda);
     129           0 :     poiss_rat = (_lambda) / ((3 * _bulk_modulus - _lambda));
     130           0 :     _effective_stiffness_local =
     131           0 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     132           0 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     133             :   }
     134       34176 :   else if (_shear_modulus_set && _youngs_modulus_set)
     135             :   {
     136       32304 :     iso_const[0] = _shear_modulus * (_youngs_modulus - 2.0 * _shear_modulus) /
     137       32304 :                    (3.0 * _shear_modulus - _youngs_modulus);
     138       32304 :     iso_const[1] = _shear_modulus;
     139       32304 :     elas_mod = _youngs_modulus;
     140       32304 :     poiss_rat = (_youngs_modulus - 2 * _shear_modulus) / (2 * _shear_modulus);
     141       32304 :     _effective_stiffness_local =
     142       64608 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     143       32304 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     144             :   }
     145        1872 :   else if (_shear_modulus_set && _poissons_ratio_set)
     146             :   {
     147        1872 :     iso_const[0] = 2.0 * _shear_modulus * _poissons_ratio / (1.0 - 2.0 * _poissons_ratio);
     148        1872 :     iso_const[1] = _shear_modulus;
     149        1872 :     elas_mod = (2 * _shear_modulus * (1 + _poissons_ratio));
     150             :     poiss_rat = (_poissons_ratio);
     151        1872 :     _effective_stiffness_local =
     152        3744 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     153        1872 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     154             :   }
     155           0 :   else if (_youngs_modulus_set && _bulk_modulus_set)
     156             :   {
     157           0 :     iso_const[0] = 3.0 * _bulk_modulus * (3.0 * _bulk_modulus - _youngs_modulus) /
     158           0 :                    (9.0 * _bulk_modulus - _youngs_modulus);
     159           0 :     iso_const[1] = 3.0 * _bulk_modulus * _youngs_modulus / (9.0 * _bulk_modulus - _youngs_modulus);
     160           0 :     elas_mod = (_youngs_modulus);
     161           0 :     poiss_rat = (3 * _bulk_modulus - _youngs_modulus) / (6 * _bulk_modulus);
     162           0 :     _effective_stiffness_local =
     163           0 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     164           0 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     165             :   }
     166           0 :   else if (_lambda_set && _poissons_ratio_set)
     167             :   {
     168           0 :     iso_const[0] = _lambda;
     169           0 :     iso_const[1] = _lambda * (1.0 - 2.0 * _poissons_ratio) / (2.0 * _poissons_ratio);
     170           0 :     elas_mod = (_lambda * (1 + _poissons_ratio) * (1 - 2 * _poissons_ratio)) / (_poissons_ratio);
     171             :     poiss_rat = (_poissons_ratio);
     172           0 :     _effective_stiffness_local =
     173           0 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     174           0 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     175             :   }
     176           0 :   else if (_lambda_set && _youngs_modulus_set)
     177             :   {
     178           0 :     iso_const[0] = _lambda;
     179           0 :     iso_const[1] = (_youngs_modulus - 3.0 * _lambda +
     180           0 :                     std::sqrt(_youngs_modulus * _youngs_modulus + 9.0 * _lambda * _lambda +
     181           0 :                               2.0 * _youngs_modulus * _lambda)) /
     182             :                    4.0;
     183           0 :     elas_mod = (_youngs_modulus);
     184           0 :     poiss_rat = (2 * _lambda) / (_youngs_modulus + _lambda +
     185           0 :                                  std::sqrt(std::pow(_youngs_modulus, 2) + 9 * std::pow(_lambda, 2) +
     186           0 :                                            2 * _youngs_modulus * _lambda));
     187           0 :     _effective_stiffness_local =
     188           0 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     189           0 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     190             :   }
     191             :   else
     192           0 :     mooseError("Incorrect combination of elastic properties in ComputeIsotropicElasticityTensor.");
     193             : 
     194             :   // Fill elasticity tensor
     195      147018 :   _Cijkl.fillFromInputVector(iso_const, T::symmetric_isotropic);
     196             : }
     197             : 
     198             : template <bool is_ad, typename T>
     199             : void
     200   171453114 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::computeQpElasticityTensor()
     201             : {
     202             :   // Assign elasticity tensor at a given quad point
     203   171453114 :   _elasticity_tensor[_qp] = _Cijkl;
     204             : 
     205             :   // Assign effective stiffness at a given quad point
     206   171453114 :   _effective_stiffness[_qp] = _effective_stiffness_local;
     207   171453114 : }
     208             : 
     209             : template class ComputeIsotropicElasticityTensorTempl<false, RankFourTensor>;
     210             : template class ComputeIsotropicElasticityTensorTempl<true, RankFourTensor>;
     211             : template class ComputeIsotropicElasticityTensorTempl<false, SymmetricRankFourTensor>;
     212             : template class ComputeIsotropicElasticityTensorTempl<true, SymmetricRankFourTensor>;