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 "PolycrystalDiffusivityTensorBase.h"
11 : #include "libmesh/quadrature.h"
12 :
13 : registerMooseObject("PhaseFieldApp", PolycrystalDiffusivityTensorBase);
14 :
15 : InputParameters
16 329 : PolycrystalDiffusivityTensorBase::validParams()
17 : {
18 329 : InputParameters params = Material::validParams();
19 329 : params.addClassDescription("Generates a diffusion tensor to distinguish between the bulk, grain "
20 : "boundaries, and surfaces");
21 658 : params.addRequiredCoupledVarWithAutoBuild(
22 : "v", "var_name_base", "op_num", "Array of coupled variables");
23 658 : params.addCoupledVar("T", "Temperature variable in Kelvin");
24 658 : params.addRequiredParam<Real>("D0", "Diffusion prefactor for vacancies in m^2/s");
25 658 : params.addRequiredParam<Real>("Em", "Vacancy migration energy in eV");
26 658 : params.addRequiredCoupledVar("c", "Vacancy phase variable");
27 658 : params.addParam<std::string>(
28 : "diffusivity_name", "D", "Name for the diffusivity material property");
29 658 : params.addParam<Real>("surfindex", 1.0, "Surface diffusion index weight");
30 658 : params.addParam<Real>("gbindex", 1.0, "Grain boundary diffusion index weight");
31 658 : params.addParam<Real>("bulkindex", 1.0, "Bulk diffusion index weight");
32 329 : return params;
33 0 : }
34 :
35 252 : PolycrystalDiffusivityTensorBase::PolycrystalDiffusivityTensorBase(
36 252 : const InputParameters & parameters)
37 : : DerivativeMaterialInterface<Material>(parameters),
38 252 : _T(coupledValue("T")),
39 252 : _c(coupledValue("c")),
40 252 : _grad_c(coupledGradient("c")),
41 252 : _c_name(coupledName("c", 0)),
42 504 : _diffusivity_name(getParam<std::string>("diffusivity_name")),
43 252 : _D(declareProperty<RealTensorValue>(_diffusivity_name)),
44 252 : _dDdc(isCoupledConstant("_c_name")
45 252 : ? nullptr
46 756 : : &declarePropertyDerivative<RealTensorValue>(_diffusivity_name, _c_name)),
47 252 : _dDdgradc(isCoupledConstant("_c_name")
48 252 : ? nullptr
49 1008 : : &declarePropertyDerivative<RankThreeTensor>(_diffusivity_name, "gradc")),
50 504 : _D0(getParam<Real>("D0")),
51 504 : _Em(getParam<Real>("Em")),
52 504 : _s_index(getParam<Real>("surfindex")),
53 504 : _gb_index(getParam<Real>("gbindex")),
54 504 : _b_index(getParam<Real>("bulkindex")),
55 252 : _kb(8.617343e-5), // Boltzmann constant in eV/K
56 252 : _op_num(coupledComponents("v")),
57 252 : _vals_name(_op_num),
58 252 : _dDdeta(_op_num),
59 504 : _dDdgradeta(_op_num)
60 : {
61 252 : if (_op_num == 0)
62 0 : mooseError("Model requires op_num > 0");
63 :
64 252 : _vals.resize(_op_num);
65 252 : _grad_vals.resize(_op_num);
66 828 : for (unsigned int i = 0; i < _op_num; ++i)
67 : {
68 576 : _vals[i] = &coupledValue("v", i);
69 576 : _grad_vals[i] = &coupledGradient("v", i);
70 1152 : _vals_name[i] = coupledName("v", i);
71 576 : if (!isCoupledConstant(_vals_name[i]))
72 576 : _dDdeta[i] = &declarePropertyDerivative<RealTensorValue>(_diffusivity_name, _vals_name[i]);
73 576 : if (!isCoupledConstant(_vals_name[i]))
74 576 : _dDdgradeta[i] =
75 1152 : &declarePropertyDerivative<RankThreeTensor>(_diffusivity_name, ("grad" + _vals_name[i]));
76 : }
77 252 : }
78 :
79 : void
80 2148553 : PolycrystalDiffusivityTensorBase::computeProperties()
81 : {
82 2148553 : RealTensorValue I(1, 0, 0, 0, 1, 0, 0, 0, 1);
83 :
84 9859565 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
85 : {
86 7711012 : Real c = _c[_qp];
87 7711012 : Real mc = 1.0 - c;
88 :
89 : // Compute grain boundary diffusivity and derivatives wrt order parameters
90 : RealTensorValue Dgb(0.0);
91 7711012 : std::vector<RealTensorValue> dDgbdeta(_op_num);
92 7711012 : std::vector<RankThreeTensor> dDgbdgradeta(_op_num);
93 :
94 23195460 : for (unsigned int i = 0; i < _op_num; ++i)
95 23351520 : for (unsigned int j = i + 1; j < _op_num; ++j)
96 : {
97 7867072 : RealGradient ngb = (*_grad_vals[i])[_qp] - (*_grad_vals[j])[_qp];
98 7867072 : if (ngb.norm() > 1.0e-10)
99 5026028 : ngb /= ngb.norm();
100 : else
101 : ngb = 0.0;
102 :
103 : RealTensorValue Tgb;
104 31468288 : for (unsigned int a = 0; a < 3; ++a)
105 70803648 : for (unsigned int b = a; b < 3; ++b)
106 : {
107 47202432 : Tgb(a, b) = I(a, b) - ngb(a) * ngb(b);
108 47202432 : Tgb(b, a) = I(b, a) - ngb(b) * ngb(a);
109 : }
110 :
111 7867072 : RankThreeTensor dTgbi, dTgbj;
112 7867072 : if (((*_grad_vals[i])[_qp] - (*_grad_vals[j])[_qp]).norm() > 1.0e-10)
113 : {
114 : Real detax = (*_grad_vals[i])[_qp](0) - (*_grad_vals[j])[_qp](0);
115 : Real detay = (*_grad_vals[i])[_qp](1) - (*_grad_vals[j])[_qp](1);
116 : Real detaz = (*_grad_vals[i])[_qp](2) - (*_grad_vals[j])[_qp](2);
117 5026028 : Real norm4 = pow(((*_grad_vals[i])[_qp] - (*_grad_vals[j])[_qp]).norm(), 4.0);
118 : // Derivatives wrt detai/dx
119 5026028 : dTgbi(0, 0, 0) = -2.0 * detax * (detay * detay + detaz * detaz) / norm4;
120 5026028 : dTgbi(1, 0, 0) = dTgbi(0, 1, 0) =
121 5026028 : (detax * detax * detay - detay * detay * detay - detay * detaz * detaz) / norm4;
122 5026028 : dTgbi(1, 1, 0) = 2.0 * detax * detay * detay / norm4;
123 5026028 : dTgbi(2, 0, 0) = dTgbi(0, 2, 0) =
124 5026028 : (detax * detax * detaz - detay * detay * detaz - detaz * detaz * detaz) / norm4;
125 5026028 : dTgbi(2, 1, 0) = dTgbi(1, 2, 0) = 2.0 * detax * detay * detaz / norm4;
126 5026028 : dTgbi(2, 2, 0) = 2.0 * detax * detaz * detaz / norm4;
127 : // Derivatives wrt detai/dy
128 5026028 : dTgbi(0, 0, 1) = 2.0 * detax * detax * detay / norm4;
129 5026028 : dTgbi(1, 0, 1) = dTgbi(0, 1, 1) =
130 5026028 : (-detax * detax * detax + detax * detay * detay - detax * detaz * detaz) / norm4;
131 5026028 : dTgbi(1, 1, 1) = -2.0 * detay * (detax * detax + detaz * detaz) / norm4;
132 5026028 : dTgbi(2, 0, 1) = dTgbi(0, 2, 1) = 2.0 * detax * detay * detaz / norm4;
133 5026028 : dTgbi(2, 1, 1) = dTgbi(1, 2, 1) =
134 5026028 : (detay * detay * detaz - detax * detax * detaz - detaz * detaz * detaz) / norm4;
135 5026028 : dTgbi(2, 2, 1) = 2.0 * detay * detaz * detaz / norm4;
136 :
137 : // Derivatives wrt detai/dz
138 5026028 : dTgbi(0, 0, 2) = 2.0 * detax * detax * detaz / norm4;
139 5026028 : dTgbi(1, 0, 2) = dTgbi(0, 1, 2) = 2.0 * detax * detay * detaz / norm4;
140 5026028 : dTgbi(1, 1, 2) = 2.0 * detay * detay * detaz / norm4;
141 5026028 : dTgbi(2, 0, 2) = dTgbi(0, 2, 2) =
142 5026028 : (detax * detaz * detaz - detax * detax * detax - detay * detay * detax) / norm4;
143 5026028 : dTgbi(2, 1, 2) = dTgbi(1, 2, 2) =
144 5026028 : (detay * detaz * detaz - detax * detax * detay - detay * detay * detay) / norm4;
145 5026028 : dTgbi(2, 2, 2) = -2.0 * detaz * (detax * detax + detay * detay) / norm4;
146 :
147 5026028 : dTgbj = -dTgbi;
148 : }
149 :
150 7867072 : Dgb += (*_vals[i])[_qp] * (*_vals[j])[_qp] * Tgb;
151 7867072 : Dgb += (*_vals[j])[_qp] * (*_vals[i])[_qp] * Tgb;
152 7867072 : dDgbdeta[i] += 2.0 * (*_vals[j])[_qp] * Tgb;
153 7867072 : dDgbdeta[j] += 2.0 * (*_vals[i])[_qp] * Tgb;
154 7867072 : dDgbdgradeta[i] += 2.0 * (*_vals[i])[_qp] * (*_vals[j])[_qp] * dTgbi;
155 7867072 : dDgbdgradeta[j] += 2.0 * (*_vals[i])[_qp] * (*_vals[j])[_qp] * dTgbj;
156 : }
157 :
158 : // Compute surface diffusivity matrix
159 : RealGradient ns(0);
160 7711012 : if (_grad_c[_qp].norm() > 1.0e-10)
161 6608370 : ns = _grad_c[_qp] / _grad_c[_qp].norm();
162 :
163 : RealTensorValue Ts;
164 30844048 : for (unsigned int a = 0; a < 3; ++a)
165 92532144 : for (unsigned int b = 0; b < 3; ++b)
166 : {
167 69399108 : Ts(a, b) = I(a, b) - ns(a) * ns(b);
168 : }
169 :
170 7711012 : RankThreeTensor dTs;
171 7711012 : if (_grad_c[_qp].norm() > 1.0e-10)
172 : {
173 6608370 : Real dcx = _grad_c[_qp](0);
174 6608370 : Real dcy = _grad_c[_qp](1);
175 6608370 : Real dcz = _grad_c[_qp](2);
176 6608370 : Real norm4 = pow(_grad_c[_qp].norm(), 4.0);
177 : // Derivatives wrt dc/dx
178 6608370 : dTs(0, 0, 0) = -2.0 * dcx * (dcy * dcy + dcz * dcz) / norm4;
179 6608370 : dTs(1, 0, 0) = dTs(0, 1, 0) = (dcx * dcx * dcy - dcy * dcy * dcy - dcy * dcz * dcz) / norm4;
180 6608370 : dTs(1, 1, 0) = 2.0 * dcx * dcy * dcy / norm4;
181 6608370 : dTs(2, 0, 0) = dTs(0, 2, 0) = (dcx * dcx * dcz - dcy * dcy * dcz - dcz * dcz * dcz) / norm4;
182 6608370 : dTs(2, 1, 0) = dTs(1, 2, 0) = 2.0 * dcx * dcy * dcz / norm4;
183 6608370 : dTs(2, 2, 0) = 2.0 * dcx * dcz * dcz / norm4;
184 :
185 : // Derivatives wrt dc/dy
186 6608370 : dTs(0, 0, 1) = 2.0 * dcx * dcx * dcy / norm4;
187 6608370 : dTs(1, 0, 1) = dTs(0, 1, 1) = (-dcx * dcx * dcx + dcx * dcy * dcy - dcx * dcz * dcz) / norm4;
188 6608370 : dTs(1, 1, 1) = -2.0 * dcy * (dcx * dcx + dcz * dcz) / norm4;
189 6608370 : dTs(2, 0, 1) = dTs(0, 2, 1) = 2.0 * dcx * dcy * dcz / norm4;
190 6608370 : dTs(2, 1, 1) = dTs(1, 2, 1) = (dcy * dcy * dcz - dcx * dcx * dcz - dcz * dcz * dcz) / norm4;
191 6608370 : dTs(2, 2, 1) = 2.0 * dcy * dcz * dcz / norm4;
192 :
193 : // Derivatives wrt dc/dz
194 6608370 : dTs(0, 0, 2) = 2.0 * dcx * dcx * dcz / norm4;
195 6608370 : dTs(1, 0, 2) = dTs(0, 1, 2) = 2.0 * dcx * dcy * dcz / norm4;
196 6608370 : dTs(1, 1, 2) = 2.0 * dcy * dcy * dcz / norm4;
197 6608370 : dTs(2, 0, 2) = dTs(0, 2, 2) = (dcx * dcz * dcz - dcx * dcx * dcx - dcy * dcy * dcx) / norm4;
198 6608370 : dTs(2, 1, 2) = dTs(1, 2, 2) = (dcy * dcz * dcz - dcx * dcx * dcy - dcy * dcy * dcy) / norm4;
199 6608370 : dTs(2, 2, 2) = -2.0 * dcz * (dcx * dcx + dcy * dcy) / norm4;
200 : }
201 :
202 7711012 : RealTensorValue Dsurf = c * c * mc * mc * Ts;
203 7711012 : RealTensorValue dDsurfdc = (2.0 * c * mc * mc - 2.0 * c * c * mc) * Ts;
204 : RankThreeTensor dDsurfdgradc = c * c * mc * mc * dTs;
205 :
206 : // Compute bulk properties
207 7711012 : _Dbulk = _D0 * std::exp(-_Em / _kb / _T[_qp]);
208 : const Real mult_bulk = 1.0;
209 : const Real dmult_bulk = 0.0;
210 :
211 : // Compute diffusion tensor
212 7711012 : _D[_qp] = _Dbulk * (_b_index * mult_bulk * I + _gb_index * Dgb + _s_index * Dsurf);
213 7711012 : if (_dDdc)
214 7711012 : (*_dDdc)[_qp] = _Dbulk * (_b_index * dmult_bulk * I + _s_index * dDsurfdc);
215 7711012 : if (_dDdgradc)
216 7711012 : (*_dDdgradc)[_qp] = _Dbulk * _s_index * dDsurfdgradc;
217 23195460 : for (unsigned int i = 0; i < _op_num; ++i)
218 : {
219 15484448 : if (_dDdeta[i])
220 15484448 : (*_dDdeta[i])[_qp] = _Dbulk * _gb_index * dDgbdeta[i];
221 15484448 : if (_dDdgradeta[i])
222 15484448 : (*_dDdgradeta[i])[_qp] = _Dbulk * _gb_index * dDgbdgradeta[i];
223 : }
224 7711012 : }
225 2148553 : }
|
