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 "PolycrystalDiffusivity.h"
      11             : 
      12             : registerMooseObject("PhaseFieldApp", PolycrystalDiffusivity);
      13             : 
      14             : InputParameters
      15          94 : PolycrystalDiffusivity::validParams()
      16             : {
      17          94 :   InputParameters params = Material::validParams();
      18          94 :   params.addClassDescription(
      19             :       "Generates a diffusion coefficient to distinguish between the bulk, pore, grain "
      20             :       "boundaries, and surfaces");
      21         188 :   params.addRequiredCoupledVarWithAutoBuild(
      22             :       "v", "var_name_base", "op_num", "Array of coupled variables");
      23         188 :   params.addParam<Real>("Dbulk", 1.0, "Diffusion coefficient for volumetric diffusion in solid");
      24         188 :   params.addParam<Real>(
      25         188 :       "Dvoid", 1.0, "Diffusion coefficient for diffusion within void/pore/bubble ");
      26         188 :   params.addParam<Real>("Dsurf", 1.0, "Diffusion coefficient for surface diffusion");
      27         188 :   params.addParam<Real>("Dgb", 1.0, "Diffusion coefficient for grain boundary diffusion");
      28         188 :   params.addRequiredCoupledVar("c", "Vacancy phase variable");
      29         188 :   params.addParam<Real>("surf_weight", 1.0, "Surface diffusion weight");
      30         188 :   params.addParam<Real>("gb_weight", 1.0, "Grain boundary diffusion weight");
      31         188 :   params.addParam<Real>("bulk_weight", 1.0, "Bulk diffusion weight");
      32         188 :   params.addParam<Real>("void_weight", 1.0, "Void diffusion weight");
      33         188 :   params.addParam<MaterialPropertyName>(
      34             :       "void_switch", "hb", "Switching Function Materials for the void");
      35         188 :   params.addParam<MaterialPropertyName>(
      36             :       "solid_switch", "hm", "Switching Function Materials for the solid");
      37         188 :   params.addParam<MaterialPropertyName>("diffusivity", "D", "The name of the diffusion property");
      38          94 :   return params;
      39           0 : }
      40             : 
      41          72 : PolycrystalDiffusivity::PolycrystalDiffusivity(const InputParameters & parameters)
      42             :   : DerivativeMaterialInterface<Material>(parameters),
      43          72 :     _c(coupledValue("c")),
      44          72 :     _c_name(coupledName("c", 0)),
      45          72 :     _op_num(coupledComponents("v")),
      46          72 :     _vals(_op_num),
      47          72 :     _diff_name(getParam<MaterialPropertyName>("diffusivity")),
      48          72 :     _diff(declareProperty<Real>(_diff_name)),
      49          72 :     _dDdc(isCoupledConstant(_c_name) ? nullptr
      50         216 :                                      : &declarePropertyDerivative<Real>(_diff_name, _c_name)),
      51         144 :     _hb(getMaterialProperty<Real>("void_switch")),
      52         144 :     _hm(getMaterialProperty<Real>("solid_switch")),
      53          72 :     _dhbdc(getMaterialPropertyDerivative<Real>("void_switch", _c_name)),
      54          72 :     _dhmdc(getMaterialPropertyDerivative<Real>("solid_switch", _c_name)),
      55          72 :     _dhbdv(_op_num),
      56          72 :     _dhmdv(_op_num),
      57         144 :     _diff_bulk(getParam<Real>("Dbulk")),
      58         144 :     _diff_void(getParam<Real>("Dvoid")),
      59         144 :     _diff_surf(getParam<Real>("Dsurf")),
      60         144 :     _diff_gb(getParam<Real>("Dgb")),
      61         144 :     _s_weight(getParam<Real>("surf_weight")),
      62         144 :     _gb_weight(getParam<Real>("gb_weight")),
      63         144 :     _b_weight(getParam<Real>("bulk_weight")),
      64         216 :     _v_weight(getParam<Real>("void_weight"))
      65             : {
      66          72 :   if (_op_num == 0)
      67           0 :     paramError("op_num", "Model requires a non zero number of order parameters.");
      68             : 
      69          72 :   _dDdv.resize(_op_num);
      70         360 :   for (MooseIndex(_op_num) op_index = 0; op_index < _op_num; ++op_index)
      71             :   {
      72         288 :     _vals[op_index] = &coupledValue("v", op_index);
      73         288 :     const VariableName op_name = coupledName("v", op_index);
      74         576 :     if (!isCoupledConstant("v"))
      75         576 :       _dDdv[op_index] = &declarePropertyDerivative<Real>(_diff_name, coupledName("v", op_index));
      76         288 :     _dhbdv[op_index] = &getMaterialPropertyDerivative<Real>("void_switch", op_index);
      77         288 :     _dhmdv[op_index] = &getMaterialPropertyDerivative<Real>("solid_switch", op_index);
      78             :   }
      79          72 : }
      80             : 
      81             : void
      82      400000 : PolycrystalDiffusivity::computeQpProperties()
      83             : {
      84             :   Real SumEtaij = 0.0;
      85             :   Real SumEtaj = 0.0;
      86     2000000 :   for (const auto i : make_range(_op_num))
      87     4000000 :     for (const auto j : make_range(i + 1, _op_num))
      88             :     {
      89     2400000 :       SumEtaij += 18.0 * (*_vals[i])[_qp] * (*_vals[i])[_qp] * (*_vals[j])[_qp] * (*_vals[j])[_qp];
      90     2400000 :       SumEtaj += 18.0 * ((*_vals[i])[_qp] * (*_vals[j])[_qp] * (*_vals[j])[_qp] +
      91     2400000 :                          (*_vals[j])[_qp] * (*_vals[i])[_qp] * (*_vals[i])[_qp]);
      92             :     }
      93      400000 :   Real c = _c[_qp];
      94      400000 :   c = std::abs(c) > 1.0 ? 1.0 : (c < 0.0 ? std::abs(c) : c);
      95      400000 :   const Real mc = 1.0 - c;
      96             : 
      97             :   // Compute diffusion function
      98      400000 :   _diff[_qp] = _b_weight * _diff_bulk * _hm[_qp] + _v_weight * _diff_void * _hb[_qp] +
      99      400000 :                30.0 * _diff_surf * _s_weight * c * c * mc * mc + _diff_gb * SumEtaij * _gb_weight;
     100             : 
     101      400000 :   if (_dDdc)
     102      400000 :     (*_dDdc)[_qp] = _b_weight * _diff_bulk * _dhmdc[_qp] + _v_weight * _diff_void * _dhbdc[_qp] +
     103      400000 :                     30.0 * _diff_surf * _s_weight * (2.0 * c * mc * mc - 2.0 * c * c * mc);
     104     2000000 :   for (const auto op_index : make_range(_op_num))
     105     1600000 :     if (_dDdv[op_index])
     106     1600000 :       (*_dDdv[op_index])[_qp] = _b_weight * _diff_bulk * (*_dhmdv[op_index])[_qp] +
     107     1600000 :                                 _v_weight * _diff_void * (*_dhbdv[op_index])[_qp] +
     108     1600000 :                                 _diff_gb * SumEtaj * _gb_weight;
     109      400000 : }