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