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 "ComputeIsotropicElasticityTensor.h"
      11             : 
      12             : registerMooseObject("TensorMechanicsApp", ComputeIsotropicElasticityTensor);
      13             : registerMooseObject("TensorMechanicsApp", ADComputeIsotropicElasticityTensor);
      14             : registerMooseObject("TensorMechanicsApp", SymmetricIsotropicElasticityTensor);
      15             : registerMooseObject("TensorMechanicsApp", ADSymmetricIsotropicElasticityTensor);
      16             : 
      17             : template <bool is_ad, typename T>
      18             : InputParameters
      19       10746 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::validParams()
      20             : {
      21       10746 :   InputParameters params = ComputeElasticityTensorBase::validParams();
      22       10746 :   params.addClassDescription("Compute a constant isotropic elasticity tensor.");
      23       21492 :   params.addParam<Real>("bulk_modulus", -1, "The bulk modulus for the material.");
      24       21492 :   params.addParam<Real>("lambda", -1, "Lame's first constant for the material.");
      25       21492 :   params.addParam<Real>("poissons_ratio", -1, "Poisson's ratio for the material.");
      26       21492 :   params.addParam<Real>("shear_modulus", -1, "The shear modulus of the material.");
      27       21492 :   params.addParam<Real>("youngs_modulus", -1, "Young's modulus of the material.");
      28       21492 :   params.declareControllable("bulk_modulus lambda poissons_ratio shear_modulus youngs_modulus");
      29       10746 :   return params;
      30           0 : }
      31             : 
      32             : template <bool is_ad, typename T>
      33        8054 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::ComputeIsotropicElasticityTensorTempl(
      34             :     const InputParameters & parameters)
      35             :   : ComputeElasticityTensorBaseTempl<is_ad, T>(parameters),
      36        8054 :     _bulk_modulus_set(parameters.isParamSetByUser("bulk_modulus")),
      37        8054 :     _lambda_set(parameters.isParamSetByUser("lambda")),
      38        8054 :     _poissons_ratio_set(parameters.isParamSetByUser("poissons_ratio")),
      39        8054 :     _shear_modulus_set(parameters.isParamSetByUser("shear_modulus")),
      40        8054 :     _youngs_modulus_set(parameters.isParamSetByUser("youngs_modulus")),
      41       16108 :     _bulk_modulus(this->template getParam<Real>("bulk_modulus")),
      42       16108 :     _lambda(this->template getParam<Real>("lambda")),
      43       16108 :     _poissons_ratio(this->template getParam<Real>("poissons_ratio")),
      44       16108 :     _shear_modulus(this->template getParam<Real>("shear_modulus")),
      45       16108 :     _youngs_modulus(this->template getParam<Real>("youngs_modulus")),
      46       16108 :     _effective_stiffness_local(parameters.isParamValid("effective_stiffness_local"))
      47             : {
      48        8054 :   unsigned int num_elastic_constants = _bulk_modulus_set + _lambda_set + _poissons_ratio_set +
      49        8054 :                                        _shear_modulus_set + _youngs_modulus_set;
      50        8054 :   if (num_elastic_constants != 2)
      51           2 :     mooseError("Exactly two elastic constants must be defined for material '" + name() + "'.");
      52             : 
      53             :   // all tensors created by this class are always isotropic
      54       16106 :   issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);
      55        8053 :   issueGuarantee("effective_stiffness", Guarantee::ISOTROPIC);
      56       16106 :   if (!isParamValid("elasticity_tensor_prefactor"))
      57       15998 :     issueGuarantee(_elasticity_tensor_name, Guarantee::CONSTANT_IN_TIME);
      58             : 
      59        8053 :   if (_bulk_modulus_set && _bulk_modulus <= 0.0)
      60           2 :     mooseError("Bulk modulus must be positive in material '" + name() + "'.");
      61             : 
      62        8052 :   if (_poissons_ratio_set && (_poissons_ratio <= -1.0 || _poissons_ratio >= 0.5))
      63           2 :     mooseError("Poissons ratio must be greater than -1 and less than 0.5 in "
      64             :                "material '" +
      65             :                name() + "'.");
      66             : 
      67        8051 :   if (_shear_modulus_set && _shear_modulus < 0.0)
      68           2 :     mooseError("Shear modulus must not be negative in material '" + name() + "'.");
      69             : 
      70        8050 :   if (_youngs_modulus_set && _youngs_modulus <= 0.0)
      71           2 :     mooseError("Youngs modulus must be positive in material '" + name() + "'.");
      72        8049 : }
      73             : 
      74             : template <bool is_ad, typename T>
      75             : void
      76     1591626 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::residualSetup()
      77             : {
      78     1591626 :   std::vector<Real> iso_const(2);
      79             :   Real elas_mod;
      80             :   Real poiss_rat;
      81             : 
      82     1591626 :   if (_youngs_modulus_set && _poissons_ratio_set)
      83             :   {
      84     1503750 :     _Cijkl.fillSymmetricIsotropicEandNu(_youngs_modulus, _poissons_ratio);
      85     1503750 :     _effective_stiffness_local =
      86     3007500 :         std::max(std::sqrt((_youngs_modulus * (1 - _poissons_ratio)) /
      87     1503750 :                            ((1 + _poissons_ratio) * (1 - 2 * _poissons_ratio))),
      88     1503750 :                  std::sqrt(_youngs_modulus / (2 * (1 + _poissons_ratio))));
      89             :     return;
      90             :   }
      91             : 
      92       87876 :   if (_lambda_set && _shear_modulus_set)
      93             :   {
      94       69639 :     iso_const[0] = _lambda;
      95       69639 :     iso_const[1] = _shear_modulus;
      96       69639 :     elas_mod = (_shear_modulus * (3 * _lambda + 2 * _shear_modulus)) / (_lambda + _shear_modulus);
      97       69639 :     poiss_rat = _lambda / (2 * (_lambda + _shear_modulus));
      98       69639 :     _effective_stiffness_local =
      99      139278 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     100       69639 :                  std::sqrt(_shear_modulus));
     101             :   }
     102       18237 :   else if (_shear_modulus_set && _bulk_modulus_set)
     103             :   {
     104         342 :     iso_const[0] = _bulk_modulus - 2.0 / 3.0 * _shear_modulus;
     105         342 :     iso_const[1] = _shear_modulus;
     106         342 :     elas_mod = (9 * _bulk_modulus * _shear_modulus) / (3 * _bulk_modulus + _shear_modulus);
     107         342 :     poiss_rat =
     108         342 :         (3 * _bulk_modulus - 2 * _shear_modulus) / (2 * (3 * _bulk_modulus + _shear_modulus));
     109         342 :     _effective_stiffness_local =
     110         684 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     111         342 :                  std::sqrt(_shear_modulus));
     112             :   }
     113       17895 :   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       17895 :   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       17895 :   else if (_shear_modulus_set && _youngs_modulus_set)
     135             :   {
     136       16167 :     iso_const[0] = _shear_modulus * (_youngs_modulus - 2.0 * _shear_modulus) /
     137       16167 :                    (3.0 * _shear_modulus - _youngs_modulus);
     138       16167 :     iso_const[1] = _shear_modulus;
     139       16167 :     elas_mod = _youngs_modulus;
     140       16167 :     poiss_rat = (_youngs_modulus - 2 * _shear_modulus) / (2 * _shear_modulus);
     141       16167 :     _effective_stiffness_local =
     142       32334 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     143       16167 :                  std::sqrt(elas_mod / (2 * (1 + poiss_rat))));
     144             :   }
     145        1728 :   else if (_shear_modulus_set && _poissons_ratio_set)
     146             :   {
     147        1728 :     iso_const[0] = 2.0 * _shear_modulus * _poissons_ratio / (1.0 - 2.0 * _poissons_ratio);
     148        1728 :     iso_const[1] = _shear_modulus;
     149        1728 :     elas_mod = (2 * _shear_modulus * (1 + _poissons_ratio));
     150             :     poiss_rat = (_poissons_ratio);
     151        1728 :     _effective_stiffness_local =
     152        3456 :         std::max(std::sqrt((elas_mod * (1 - poiss_rat)) / ((1 + poiss_rat) * (1 - 2 * poiss_rat))),
     153        1728 :                  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       87876 :   _Cijkl.fillFromInputVector(iso_const, T::symmetric_isotropic);
     196             : }
     197             : 
     198             : template <bool is_ad, typename T>
     199             : void
     200    84423019 : ComputeIsotropicElasticityTensorTempl<is_ad, T>::computeQpElasticityTensor()
     201             : {
     202             :   // Assign elasticity tensor at a given quad point
     203    84423019 :   _elasticity_tensor[_qp] = _Cijkl;
     204             : 
     205             :   // Assign effective stiffness at a given quad point
     206    84423019 :   _effective_stiffness[_qp] = _effective_stiffness_local;
     207    84423019 : }
     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>;