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 "GBAnisotropyBase.h"
      11             : #include "MooseMesh.h"
      12             : 
      13             : #include <fstream>
      14             : 
      15             : InputParameters
      16         188 : GBAnisotropyBase::validParams()
      17             : {
      18         188 :   InputParameters params = Material::validParams();
      19         376 :   params.addCoupledVar("T", 300.0, "Temperature in Kelvin");
      20         376 :   params.addParam<Real>("length_scale", 1.0e-9, "Length scale in m, where default is nm");
      21         376 :   params.addParam<Real>("time_scale", 1.0e-9, "Time scale in s, where default is ns");
      22         376 :   params.addParam<Real>("molar_volume_value",
      23         376 :                         7.11e-6,
      24             :                         "molar volume of material in m^3/mol, by default it's the value of copper");
      25         376 :   params.addParam<Real>(
      26         376 :       "delta_sigma", 0.1, "factor determining inclination dependence of GB energy");
      27         376 :   params.addParam<Real>(
      28         376 :       "delta_mob", 0.1, "factor determining inclination dependence of GB mobility");
      29         376 :   params.addRequiredParam<FileName>("Anisotropic_GB_file_name",
      30             :                                     "Name of the file containing: 1)GB mobility prefactor; 2) GB "
      31             :                                     "migration activation energy; 3)GB energy");
      32         376 :   params.addRequiredParam<bool>("inclination_anisotropy",
      33             :                                 "The GB anisotropy inclination would be considered if true");
      34         376 :   params.addRequiredCoupledVarWithAutoBuild(
      35             :       "v", "var_name_base", "op_num", "Array of coupled variables");
      36         188 :   return params;
      37           0 : }
      38             : 
      39         144 : GBAnisotropyBase::GBAnisotropyBase(const InputParameters & parameters)
      40             :   : Material(parameters),
      41         288 :     _mesh_dimension(_mesh.dimension()),
      42         288 :     _length_scale(getParam<Real>("length_scale")),
      43         288 :     _time_scale(getParam<Real>("time_scale")),
      44         288 :     _M_V(getParam<Real>("molar_volume_value")),
      45         288 :     _delta_sigma(getParam<Real>("delta_sigma")),
      46         288 :     _delta_mob(getParam<Real>("delta_mob")),
      47         144 :     _Anisotropic_GB_file_name(getParam<FileName>("Anisotropic_GB_file_name")),
      48         288 :     _inclination_anisotropy(getParam<bool>("inclination_anisotropy")),
      49         144 :     _T(coupledValue("T")),
      50         144 :     _kappa(declareProperty<Real>("kappa_op")),
      51         144 :     _gamma(declareProperty<Real>("gamma_asymm")),
      52         144 :     _L(declareProperty<Real>("L")),
      53         144 :     _mu(declareProperty<Real>("mu")),
      54         144 :     _molar_volume(declareProperty<Real>("molar_volume")),
      55         144 :     _entropy_diff(declareProperty<Real>("entropy_diff")),
      56         144 :     _act_wGB(declareProperty<Real>("act_wGB")),
      57         144 :     _kb(8.617343e-5),      // Boltzmann constant in eV/K
      58         144 :     _JtoeV(6.24150974e18), // Joule to eV conversion
      59         144 :     _mu_qp(0.0),
      60         144 :     _op_num(coupledComponents("v")),
      61         144 :     _vals(coupledValues("v")),
      62         288 :     _grad_vals(coupledGradients("v"))
      63             : {
      64             :   // reshape vectors
      65         144 :   _sigma.resize(_op_num);
      66         144 :   _mob.resize(_op_num);
      67         144 :   _Q.resize(_op_num);
      68         144 :   _kappa_gamma.resize(_op_num);
      69         144 :   _a_g2.resize(_op_num);
      70             : 
      71         540 :   for (unsigned int op = 0; op < _op_num; ++op)
      72             :   {
      73         396 :     _sigma[op].resize(_op_num);
      74         396 :     _mob[op].resize(_op_num);
      75         396 :     _Q[op].resize(_op_num);
      76         396 :     _kappa_gamma[op].resize(_op_num);
      77         396 :     _a_g2[op].resize(_op_num);
      78             :   }
      79             : 
      80             :   // Read in data from "Anisotropic_GB_file_name"
      81         144 :   std::ifstream inFile(_Anisotropic_GB_file_name.c_str());
      82             : 
      83         144 :   if (!inFile)
      84           0 :     paramError("Anisotropic_GB_file_name", "Can't open GB anisotropy input file");
      85             : 
      86         432 :   for (unsigned int i = 0; i < 2; ++i)
      87         288 :     inFile.ignore(255, '\n'); // ignore line
      88             : 
      89             :   Real data;
      90        1332 :   for (unsigned int i = 0; i < 3 * _op_num; ++i)
      91             :   {
      92             :     std::vector<Real> row; // create an empty row of double values
      93        4536 :     for (unsigned int j = 0; j < _op_num; ++j)
      94             :     {
      95             :       inFile >> data;
      96        3348 :       row.push_back(data);
      97             :     }
      98             : 
      99        1188 :     if (i < _op_num)
     100         396 :       _sigma[i] = row; // unit: J/m^2
     101             : 
     102         792 :     else if (i < 2 * _op_num)
     103         396 :       _mob[i - _op_num] = row; // unit: m^4/(J*s)
     104             : 
     105             :     else
     106         396 :       _Q[i - 2 * _op_num] = row; // unit: eV
     107        1188 :   }
     108             : 
     109         144 :   inFile.close();
     110         144 : }
     111             : 
     112             : void
     113     5925360 : GBAnisotropyBase::computeQpProperties()
     114             : {
     115             :   Real sum_kappa = 0.0;
     116             :   Real sum_gamma = 0.0;
     117             :   Real sum_L = 0.0;
     118             :   Real Val = 0.0;
     119             :   Real sum_val = 0.0;
     120             :   Real f_sigma = 1.0;
     121             :   Real f_mob = 1.0;
     122             :   Real gamma_value = 0.0;
     123             : 
     124    15491280 :   for (unsigned int m = 0; m < _op_num - 1; ++m)
     125             :   {
     126    22772400 :     for (unsigned int n = m + 1; n < _op_num; ++n) // m<n
     127             :     {
     128    13206480 :       gamma_value = _kappa_gamma[n][m];
     129             : 
     130    13206480 :       if (_inclination_anisotropy)
     131             :       {
     132     2284800 :         if (_mesh_dimension == 3)
     133           0 :           mooseError("This material doesn't support inclination dependence for 3D for now!");
     134             : 
     135     2284800 :         Real phi_ave = libMesh::pi * n / (2.0 * _op_num);
     136     2284800 :         Real sin_phi = std::sin(2.0 * phi_ave);
     137     2284800 :         Real cos_phi = std::cos(2.0 * phi_ave);
     138             : 
     139     2284800 :         Real a = (*_grad_vals[m])[_qp](0) - (*_grad_vals[n])[_qp](0);
     140     2284800 :         Real b = (*_grad_vals[m])[_qp](1) - (*_grad_vals[n])[_qp](1);
     141     2284800 :         Real ab = a * a + b * b + 1.0e-7; // for the sake of numerical convergence, the smaller the
     142             :                                           // more accurate, but more difficult to converge
     143             : 
     144     2284800 :         Real cos_2phi = cos_phi * (a * a - b * b) / ab + sin_phi * 2.0 * a * b / ab;
     145     2284800 :         Real cos_4phi = 2.0 * cos_2phi * cos_2phi - 1.0;
     146             : 
     147     2284800 :         f_sigma = 1.0 + _delta_sigma * cos_4phi;
     148     2284800 :         f_mob = 1.0 + _delta_mob * cos_4phi;
     149             : 
     150     2284800 :         Real g2 = _a_g2[n][m] * f_sigma;
     151     2284800 :         Real y = -5.288 * g2 * g2 * g2 * g2 - 0.09364 * g2 * g2 * g2 + 9.965 * g2 * g2 -
     152     2284800 :                  8.183 * g2 + 2.007;
     153     2284800 :         gamma_value = 1.0 / y;
     154             :       }
     155             : 
     156    13206480 :       Val = (100000.0 * ((*_vals[m])[_qp]) * ((*_vals[m])[_qp]) + 0.01) *
     157    13206480 :             (100000.0 * ((*_vals[n])[_qp]) * ((*_vals[n])[_qp]) + 0.01);
     158             : 
     159    13206480 :       sum_val += Val;
     160    13206480 :       sum_kappa += _kappa_gamma[m][n] * f_sigma * Val;
     161    13206480 :       sum_gamma += gamma_value * Val;
     162             :       // Following comes from substituting Eq. (36c) from the paper into (36b)
     163    13206480 :       sum_L += Val * _mob[m][n] * std::exp(-_Q[m][n] / (_kb * _T[_qp])) * f_mob * _mu_qp *
     164    13206480 :                _a_g2[n][m] / _sigma[m][n];
     165             :     }
     166             :   }
     167             : 
     168     5925360 :   _kappa[_qp] = sum_kappa / sum_val;
     169     5925360 :   _gamma[_qp] = sum_gamma / sum_val;
     170     5925360 :   _L[_qp] = sum_L / sum_val;
     171     5925360 :   _mu[_qp] = _mu_qp;
     172             : 
     173     5925360 :   _molar_volume[_qp] =
     174     5925360 :       _M_V / (_length_scale * _length_scale * _length_scale); // m^3/mol converted to ls^3/mol
     175     5925360 :   _entropy_diff[_qp] = 9.5 * _JtoeV;                          // J/(K mol) converted to eV(K mol)
     176     5925360 :   _act_wGB[_qp] = 0.5e-9 / _length_scale;                     // 0.5 nm
     177     5925360 : }