Line data    Source code

       1             : /*************************************************/
       2             : /*           DO NOT MODIFY THIS HEADER           */
       3             : /*                                               */
       4             : /*                     MASTODON                  */
       5             : /*                                               */
       6             : /*    (c) 2015 Battelle Energy Alliance, LLC     */
       7             : /*            ALL RIGHTS RESERVED                */
       8             : /*                                               */
       9             : /*   Prepared by Battelle Energy Alliance, LLC   */
      10             : /*     With the U. S. Department of Energy       */
      11             : /*                                               */
      12             : /*     See COPYRIGHT for full restrictions       */
      13             : /*************************************************/
      14             : #include "ComputeIsotropicElasticityTensorSoil.h"
      15             : 
      16             : registerMooseObject("MastodonApp", ComputeIsotropicElasticityTensorSoil);
      17             : registerMooseObject("MastodonApp", ADComputeIsotropicElasticityTensorSoil);
      18             : 
      19             : template <bool is_ad>
      20             : InputParameters
      21         994 : ComputeIsotropicElasticityTensorSoilTempl<is_ad>::validParams()
      22             : {
      23         994 :   InputParameters params = ComputeElasticityTensorBaseTempl<is_ad>::validParams();
      24         994 :   params += LayeredMaterialInterface<ComputeElasticityTensorBaseTempl<is_ad>>::validParams();
      25         994 :   params.addClassDescription("Compute an isotropic elasticity tensor for a "
      26             :                              "layered soil material when shear modulus or elastic modulus, "
      27             :                              "poisson's ratio and density are provided as "
      28             :                              "input for each layer.");
      29        1988 :   params.setDocString("layer_ids",
      30             :                       "Vector of layer ids that map one-to-one "
      31             :                       "with the 'shear_modulus' or 'elastic_modulus', 'poissons_ratio' "
      32             :                       "and 'density' input parameters.");
      33        1988 :   params.addRequiredParam<std::vector<Real>>("poissons_ratio",
      34             :                                              "Vector of Poisson's ratio values that map one-to-one "
      35             :                                              "with the number 'layer_ids' parameter.");
      36        1988 :   params.addParam<std::vector<Real>>("shear_modulus",
      37             :                                      "Vector of shear modulus values that map one-to-one "
      38             :                                      "with the number 'layer_ids' parameter.");
      39        1988 :   params.addParam<std::vector<Real>>("elastic_modulus",
      40             :                                      "Vector of elastic modulus values that map one-to-one "
      41             :                                      "with the number 'layer_ids' parameter.");
      42        1988 :   params.addRequiredParam<std::vector<Real>>(
      43             :       "density",
      44             :       "Vector of density values that map one-to-one with the number "
      45             :       "'layer_ids' parameter.");
      46        1988 :   params.addParam<bool>(
      47        1988 :       "wave_speed_calculation", true, "Set to False to turn off P and S wave speed calculation.");
      48             :   // Controlled scale parameters (temporary; TODO: remove when #107 is resolved.)
      49        1988 :   params.addParam<Real>("scale_factor_density", 1.0, "Scale factor for density.");
      50        1988 :   params.declareControllable("scale_factor_density");
      51         994 :   return params;
      52           0 : }
      53             : 
      54             : template <bool is_ad>
      55             : const MooseArray<Real> &
      56         766 : getShearModulus(ComputeIsotropicElasticityTensorSoilTempl<is_ad> * object,
      57             :                 std::vector<Real> & shear_modulus,
      58             :                 std::string name)
      59             : {
      60        2470 :   if (object->isParamValid("shear_modulus") && object->isParamValid("elastic_modulus"))
      61           0 :     mooseError("In block " + name +
      62             :                ". Please provide ONE of the parameters, 'shear_modulus' and "
      63             :                "'elastic_modulus', but not both.");
      64        2423 :   if (!object->isParamValid("shear_modulus") && !object->isParamValid("elastic_modulus"))
      65           0 :     mooseError("In block " + name +
      66             :                ". Please provide ONE of the parameters, 'shear_modulus' or 'elastic_modulus'.");
      67        1532 :   if (object->isParamValid("shear_modulus"))
      68         937 :     return object->template getLayerParam<Real>("shear_modulus");
      69             :   else
      70             :   {
      71             :     const std::vector<Real> & elastic_modulus =
      72         297 :         object->template getParam<std::vector<Real>>("elastic_modulus");
      73             :     const std::vector<Real> & poissons_ratio =
      74         594 :         object->template getParam<std::vector<Real>>("poissons_ratio");
      75         297 :     shear_modulus.resize(elastic_modulus.size());
      76        1134 :     for (std::size_t i = 0; i < shear_modulus.size(); ++i)
      77         837 :       shear_modulus[i] = elastic_modulus[i] / (2 * (1 + poissons_ratio[i]));
      78         297 :     return object->addLayerVector(shear_modulus);
      79             :   }
      80             : }
      81             : 
      82             : template <bool is_ad>
      83         766 : ComputeIsotropicElasticityTensorSoilTempl<is_ad>::ComputeIsotropicElasticityTensorSoilTempl(
      84             :     const InputParameters & parameters)
      85             :   : LayeredMaterialInterface<ComputeElasticityTensorBaseTempl<is_ad>>(parameters),
      86             :     _input_shear_modulus(),
      87         766 :     _layer_shear_modulus(getShearModulus(this, _input_shear_modulus, name())),
      88         765 :     _layer_density(this->template getLayerParam<Real>("density")),
      89         765 :     _layer_poissons_ratio(this->template getLayerParam<Real>("poissons_ratio")),
      90        1530 :     _wave_speed_calculation(this->template getParam<bool>("wave_speed_calculation")),
      91        1530 :     _shear_wave_speed(_wave_speed_calculation
      92        1152 :                           ? &this->template declareGenericProperty<Real, is_ad>("shear_wave_speed")
      93             :                           : NULL),
      94        1530 :     _P_wave_speed(_wave_speed_calculation
      95        1152 :                       ? &this->template declareGenericProperty<Real, is_ad>("P_wave_speed")
      96             :                       : NULL),
      97        1530 :     _density(this->template declareGenericProperty<Real, is_ad>("density")),
      98        1530 :     _scale_density(this->template getParam<Real>("scale_factor_density")),
      99        1531 :     _effective_stiffness_local(parameters.isParamValid("effective_stiffness_local"))
     100             : {
     101             : 
     102             :   // all tensors created by this class are always isotropic
     103        1530 :   issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);
     104        1530 :   issueGuarantee("effective_stiffness", Guarantee::ISOTROPIC);
     105             : 
     106             :   // all tensors created by this class are always constant in time
     107         765 :   issueGuarantee(_elasticity_tensor_name, Guarantee::CONSTANT_IN_TIME);
     108             : 
     109         765 :   std::vector<Real> iso_const(2);
     110             :   // Fill elasticity tensor
     111         765 :   _Cijkl.fillFromInputVector(iso_const, RankFourTensor::symmetric_isotropic);
     112         765 : }
     113             : 
     114             : template <bool is_ad>
     115             : void
     116     9220552 : ComputeIsotropicElasticityTensorSoilTempl<is_ad>::computeQpElasticityTensor()
     117             : {
     118             : 
     119     9220552 :   _P_wave_modulus = _layer_shear_modulus[_qp] * 2.0 * (1.0 - _layer_poissons_ratio[_qp]) /
     120     9220552 :                     (1.0 - 2.0 * _layer_poissons_ratio[_qp]);
     121             : 
     122     9220552 :   _density[_qp] = _layer_density[_qp] * _scale_density;
     123             : 
     124     9220552 :   if (_wave_speed_calculation)
     125             :   {
     126             :     // Shear wave speed: sqrt(G/rho)
     127     6387304 :     (*_shear_wave_speed)[_qp] = std::sqrt(_layer_shear_modulus[_qp] / _density[_qp]);
     128             : 
     129             :     // P wave speed: sqrt(M/rho)
     130     6387304 :     (*_P_wave_speed)[_qp] = std::sqrt(_P_wave_modulus / _density[_qp]);
     131             :   }
     132             : 
     133     9220552 :   std::vector<Real> iso_const(2);
     134     9220552 :   iso_const[0] = _P_wave_modulus - 2.0 * _layer_shear_modulus[_qp]; // lambda = M - 2G
     135     9220552 :   iso_const[1] = _layer_shear_modulus[_qp];                         // shear modulus
     136             : 
     137             :   // Fill elasticity tensor
     138     9220552 :   _Cijkl.fillFromInputVector(iso_const, RankFourTensor::symmetric_isotropic);
     139             : 
     140             :   // Effective stiffness computations
     141     9220552 :   Real elas_mod = _layer_shear_modulus[_qp] * 2.0 * (1.0 + _layer_poissons_ratio[_qp]);
     142     9220552 :   _effective_stiffness_local =
     143    18441104 :       std::max(std::sqrt((elas_mod * (1 - _layer_poissons_ratio[_qp])) /
     144     9220552 :                          ((1 + _layer_poissons_ratio[_qp]) * (1 - 2 * _layer_poissons_ratio[_qp]))),
     145     9220552 :                std::sqrt(elas_mod / (2 * (1 + _layer_poissons_ratio[_qp]))));
     146             : 
     147             :   // Assign elasticity tensor at a given quad point
     148     9220552 :   _elasticity_tensor[_qp] = _Cijkl;
     149             : 
     150             :   // Assign effective stiffness at a given quad point
     151     9220552 :   _effective_stiffness[_qp] = _effective_stiffness_local;
     152     9220552 : }
     153             : 
     154             : template class ComputeIsotropicElasticityTensorSoilTempl<false>;
     155             : template class ComputeIsotropicElasticityTensorSoilTempl<true>;