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 "ComputePolycrystalElasticityTensor.h"
      11             : #include "RotationTensor.h"
      12             : 
      13             : registerMooseObject("PhaseFieldApp", ComputePolycrystalElasticityTensor);
      14             : 
      15             : InputParameters
      16          31 : ComputePolycrystalElasticityTensor::validParams()
      17             : {
      18          31 :   InputParameters params = ComputeElasticityTensorBase::validParams();
      19          31 :   params.addClassDescription(
      20             :       "Compute an evolving elasticity tensor coupled to a grain growth phase field model.");
      21          62 :   params.addRequiredParam<UserObjectName>(
      22             :       "grain_tracker", "Name of GrainTracker user object that provides RankFourTensors");
      23          62 :   params.addParam<Real>("length_scale", 1.0e-9, "Length scale of the problem, in meters");
      24          62 :   params.addParam<Real>("pressure_scale", 1.0e6, "Pressure scale of the problem, in pa");
      25          62 :   params.addRequiredCoupledVarWithAutoBuild(
      26             :       "v", "var_name_base", "op_num", "Array of coupled variables");
      27          62 :   params.addParam<UserObjectName>("euler_angle_provider",
      28             :                                   "Name of Euler angle provider user object");
      29             : 
      30          31 :   return params;
      31           0 : }
      32             : 
      33          24 : ComputePolycrystalElasticityTensor::ComputePolycrystalElasticityTensor(
      34          24 :     const InputParameters & parameters)
      35             :   : ComputeElasticityTensorBase(parameters),
      36          24 :     _length_scale(getParam<Real>("length_scale")),
      37          48 :     _pressure_scale(getParam<Real>("pressure_scale")),
      38          24 :     _grain_tracker(getUserObject<GrainDataTracker<RankFourTensor>>("grain_tracker")),
      39          24 :     _op_num(coupledComponents("v")),
      40          24 :     _vals(coupledValues("v")),
      41          24 :     _euler(isParamValid("euler_angle_provider")
      42          48 :                ? &getUserObject<EulerAngleProvider>("euler_angle_provider")
      43             :                : nullptr),
      44          24 :     _crysrot(isParamValid("euler_angle_provider")
      45          48 :                  ? &declareProperty<RankTwoTensor>(_base_name + "crysrot")
      46             :                  : nullptr),
      47          24 :     _D_elastic_tensor(_op_num),
      48          24 :     _JtoeV(6.24150974e18)
      49             : {
      50             :   // Loop over variables (ops)
      51         120 :   for (MooseIndex(_op_num) op_index = 0; op_index < _op_num; ++op_index)
      52             :   {
      53             :     // declare elasticity tensor derivative properties
      54         288 :     _D_elastic_tensor[op_index] = &declarePropertyDerivative<RankFourTensor>(
      55         192 :         _elasticity_tensor_name, coupledName("v", op_index));
      56             :   }
      57          24 : }
      58             : 
      59             : void
      60     1298304 : ComputePolycrystalElasticityTensor::computeQpElasticityTensor()
      61             : {
      62             :   // Get list of active order parameters from grain tracker
      63     1298304 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
      64             : 
      65             :   // Calculate elasticity tensor
      66     1298304 :   _elasticity_tensor[_qp].zero();
      67     1298304 :   Real sum_h = 0.0;
      68     6491520 :   for (MooseIndex(op_to_grains) op_index = 0; op_index < op_to_grains.size(); ++op_index)
      69             :   {
      70     5193216 :     auto grain_id = op_to_grains[op_index];
      71     5193216 :     if (grain_id == FeatureFloodCount::invalid_id)
      72             :       continue;
      73             : 
      74             :     // Interpolation factor for elasticity tensors
      75     1689120 :     Real h = (1.0 + std::sin(libMesh::pi * ((*_vals[op_index])[_qp] - 0.5))) / 2.0;
      76             : 
      77             :     // Sum all rotated elasticity tensors
      78     3378240 :     _elasticity_tensor[_qp] += _grain_tracker.getData(grain_id) * h;
      79     1689120 :     sum_h += h;
      80             : 
      81     3378240 :     if (isParamValid("euler_angle_provider"))
      82             :     {
      83     1689120 :       EulerAngles angles;
      84     1689120 :       angles = _euler->getEulerAngles(grain_id);
      85             : 
      86     1689120 :       RotationTensor R(angles);
      87     1689120 :       if ((*_vals[op_index])[_qp] > 0.5)
      88     1283532 :         (*_crysrot)[_qp] = R; // this is done for the crystal plasticity model compatibility
      89             :     }
      90             :   }
      91             : 
      92     1298304 :   const Real tol = 1.0e-10;
      93     1298304 :   sum_h = std::max(sum_h, tol);
      94     1298304 :   _elasticity_tensor[_qp] /= sum_h;
      95             : 
      96             :   // Calculate elasticity tensor derivative: Cderiv = dhdopi/sum_h * (Cop - _Cijkl)
      97     6491520 :   for (MooseIndex(_op_num) op_index = 0; op_index < _op_num; ++op_index)
      98     5193216 :     (*_D_elastic_tensor[op_index])[_qp].zero();
      99             : 
     100     6491520 :   for (MooseIndex(op_to_grains) op_index = 0; op_index < op_to_grains.size(); ++op_index)
     101             :   {
     102     5193216 :     auto grain_id = op_to_grains[op_index];
     103     5193216 :     if (grain_id == FeatureFloodCount::invalid_id)
     104             :       continue;
     105             : 
     106     1689120 :     Real dhdopi = libMesh::pi * std::cos(libMesh::pi * ((*_vals[op_index])[_qp] - 0.5)) / 2.0;
     107     1689120 :     RankFourTensor & C_deriv = (*_D_elastic_tensor[op_index])[_qp];
     108             : 
     109     1689120 :     C_deriv = (_grain_tracker.getData(grain_id) - _elasticity_tensor[_qp]) * dhdopi / sum_h;
     110             : 
     111             :     // Convert from XPa to eV/(xm)^3, where X is pressure scale and x is length scale;
     112     1689120 :     C_deriv *= _JtoeV * (_length_scale * _length_scale * _length_scale) * _pressure_scale;
     113             :   }
     114     1298304 : }