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 "ElectrochemicalSinteringMaterial.h"
11 : #include "libmesh/quadrature.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("PhaseFieldApp", ElectrochemicalSinteringMaterial);
15 :
16 : InputParameters
17 47 : ElectrochemicalSinteringMaterial::validParams()
18 : {
19 47 : InputParameters params = Material::validParams();
20 47 : params.addClassDescription("Includes switching and thermodynamic properties for the grand "
21 : "potential sintering model coupled with electrochemistry");
22 94 : params.addRequiredCoupledVarWithAutoBuild(
23 : "etas", "var_name_base", "op_num", "Array of order parameters that describe solid phase");
24 94 : params.addRequiredCoupledVar("chemical_potentials",
25 : "The name of the chemical potential variables for defects");
26 94 : params.addRequiredCoupledVar("void_op", "The name of the void phase order parameter");
27 94 : params.addRequiredCoupledVar("Temperature", "Name of the temperature variable with units of K");
28 94 : params.addRequiredCoupledVar("electric_potential",
29 : "Name of the electric potential variable with units of V");
30 94 : params.addParam<std::vector<MaterialPropertyName>>(
31 : "solid_energy_coefficients",
32 : "Vector of parabolic solid energy coefficients (energy*volume) for "
33 : "each defect species. Only used for parabolic energy. Place in same order as "
34 : "chemical_potentials!");
35 94 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
36 : "void_energy_coefficients",
37 : "Vector of parabolic void energy coefficients (energy*volume) for each defect species. Place "
38 : "in same order as chemical_potentials!");
39 94 : params.addParam<Real>("surface_energy", 19.7, "Surface energy in units of problem (energy/area)");
40 94 : params.addParam<Real>(
41 94 : "grainboundary_energy", 9.86, "Grain boundary energy in units of problem (energy/area)");
42 94 : params.addParam<Real>("int_width", 1, "Interface width in units of problem (length)");
43 141 : params.addRangeCheckedParam<Real>(
44 : "surface_switch_value",
45 94 : 0.3,
46 : "surface_switch_value >= 0 & surface_switch_value <= 1.0",
47 : "Value between 0 and 1 that determines when the interface begins to switch "
48 : "from surface to GB. Small values give less error while large values "
49 : "converge better.");
50 94 : params.addRequiredParam<std::vector<MaterialPropertyName>>(
51 : "min_vacancy_concentrations_solid",
52 : "Vector of names of materials that determine the minimum in energy wrt defect concentrations "
53 : "in the solid phase. Place in same order as chemical_potentials!");
54 94 : params.addRequiredParam<std::vector<Real>>("min_vacancy_concentrations_void",
55 : "Vector of minima in energy wrt defect concentrations "
56 : "in the void phase. Place in same order as "
57 : "chemical_potentials!");
58 94 : MooseEnum solid_energy_model("PARABOLIC DILUTE IDEAL", "PARABOLIC");
59 94 : params.addParam<MooseEnum>("solid_energy_model",
60 : solid_energy_model,
61 : "Type of energy function to use for the solid phase.");
62 94 : params.addRequiredParam<std::vector<int>>(
63 : "defect_charges",
64 : "Vector of charges of defect species. Place in same order as chemical_potentials!");
65 94 : params.addRequiredParam<Real>("solid_relative_permittivity",
66 : "Solid phase relative permittivity (dimensionless)");
67 94 : params.addParam<Real>("voltage_scale", 1, "Voltage scale (default is for voltage in V)");
68 47 : return params;
69 47 : }
70 :
71 36 : ElectrochemicalSinteringMaterial::ElectrochemicalSinteringMaterial(
72 36 : const InputParameters & parameters)
73 : : DerivativeMaterialInterface<Material>(parameters),
74 36 : _neta(coupledComponents("etas")),
75 36 : _eta(_neta),
76 36 : _eta_name(_neta),
77 36 : _ndefects(coupledComponents("chemical_potentials")),
78 36 : _w(coupledValues("chemical_potentials")),
79 36 : _w_name(_ndefects),
80 36 : _phi(coupledValue("void_op")),
81 36 : _phi_name(coupledName("void_op", 0)),
82 72 : _ns_min_names(getParam<std::vector<MaterialPropertyName>>("min_vacancy_concentrations_solid")),
83 36 : _ns_min(_ndefects),
84 36 : _dns_min(_ndefects),
85 36 : _d2ns_min(_ndefects),
86 36 : _temp(coupledValue("Temperature")),
87 36 : _v(coupledValue("electric_potential")),
88 36 : _grad_V(coupledGradient("electric_potential")),
89 72 : _kv_names(getParam<std::vector<MaterialPropertyName>>("void_energy_coefficients")),
90 36 : _kv(_ndefects),
91 72 : _ks_names(getParam<std::vector<MaterialPropertyName>>("solid_energy_coefficients")),
92 36 : _ks(_ndefects),
93 36 : _hv(declareProperty<Real>("hv")),
94 36 : _dhv(declarePropertyDerivative<Real>("hv", _phi_name)),
95 36 : _d2hv(declarePropertyDerivative<Real>("hv", _phi_name, _phi_name)),
96 36 : _hs(declareProperty<Real>("hs")),
97 36 : _dhs(declarePropertyDerivative<Real>("hs", _phi_name)),
98 36 : _d2hs(declarePropertyDerivative<Real>("hs", _phi_name, _phi_name)),
99 36 : _omegav(declareProperty<Real>("omegav")),
100 36 : _domegavdw(_ndefects),
101 36 : _d2omegavdw2(_ndefects),
102 36 : _omegas(declareProperty<Real>("omegas")),
103 36 : _domegasdw(_ndefects),
104 36 : _d2omegasdw2(_ndefects),
105 36 : _domegasdeta(_neta),
106 36 : _d2omegasdwdeta(_ndefects),
107 36 : _d2omegasdetadeta(_neta),
108 36 : _mu(declareProperty<Real>("mu")),
109 36 : _dmu(declarePropertyDerivative<Real>("mu", _phi_name)),
110 36 : _d2mu(declarePropertyDerivative<Real>("mu", _phi_name, _phi_name)),
111 36 : _kappa(declareProperty<Real>("kappa")),
112 36 : _dkappa(declarePropertyDerivative<Real>("kappa", _phi_name)),
113 36 : _d2kappa(declarePropertyDerivative<Real>("kappa", _phi_name, _phi_name)),
114 36 : _gamma(declareProperty<Real>("gamma")),
115 72 : _sigma_s(getParam<Real>("surface_energy")),
116 72 : _sigma_gb(getParam<Real>("grainboundary_energy")),
117 72 : _int_width(getParam<Real>("int_width")),
118 72 : _switch(getParam<Real>("surface_switch_value")),
119 72 : _solid_energy(getParam<MooseEnum>("solid_energy_model")),
120 36 : _mu_s(6.0 * _sigma_s / _int_width),
121 36 : _mu_gb(6.0 * _sigma_gb / _int_width),
122 36 : _kappa_s(0.75 * _sigma_s * _int_width),
123 36 : _kappa_gb(0.75 * _sigma_gb * _int_width),
124 36 : _kB(8.617343e-5), // eV/K
125 36 : _e(1.0), // to put energy units in eV
126 72 : _z(getParam<std::vector<int>>("defect_charges")),
127 72 : _v_scale(getParam<Real>("voltage_scale")),
128 72 : _eps_r(getParam<Real>("solid_relative_permittivity")),
129 108 : _nv_min(getParam<std::vector<Real>>("min_vacancy_concentrations_void"))
130 : {
131 36 : if (_ndefects != _z.size())
132 0 : paramError("defect_charges", "Need to pass in as many defect_charges as chemical_potentials");
133 36 : if (_ndefects != _ns_min_names.size())
134 0 : paramError("min_vacancy_concentrations_solid",
135 : "Need to pass in as many min_vacancy_concentrations_solid as chemical_potentials");
136 36 : if (_ndefects != _nv_min.size())
137 0 : paramError("min_vacancy_concentrations_void",
138 : "Need to pass in as many min_vacancy_concentrations_void as chemical_potentials");
139 36 : if (_ndefects != _kv_names.size())
140 0 : paramError("void_energy_coefficients",
141 : "Need to pass in as many void_energy_coefficients as chemical_potentials");
142 36 : if (_ndefects != _ks_names.size())
143 0 : paramError("solid_energy_coefficients",
144 : "Need to pass in as many solid_energy_coefficients as chemical_potentials");
145 :
146 108 : for (unsigned int i = 0; i < _ndefects; ++i)
147 : {
148 72 : _w_name[i] = coupledName("chemical_potentials", i);
149 72 : _ns_min[i] = &getMaterialPropertyByName<Real>(_ns_min_names[i]);
150 72 : _dns_min[i].resize(_neta);
151 72 : _d2ns_min[i].resize(_neta);
152 72 : _d2omegasdwdeta[i].resize(_neta);
153 72 : _kv[i] = &getMaterialPropertyByName<Real>(_kv_names[i]);
154 72 : _ks[i] = &getMaterialPropertyByName<Real>(_ks_names[i]);
155 72 : _domegavdw[i] = &declarePropertyDerivative<Real>("omegav", _w_name[i]);
156 72 : _d2omegavdw2[i] = &declarePropertyDerivative<Real>("omegav", _w_name[i], _w_name[i]);
157 72 : _domegasdw[i] = &declarePropertyDerivative<Real>("omegas", _w_name[i]);
158 72 : _d2omegasdw2[i] = &declarePropertyDerivative<Real>("omegas", _w_name[i], _w_name[i]);
159 216 : for (unsigned int j = 0; j < _neta; ++j)
160 : {
161 144 : _dns_min[i][j] = &getMaterialPropertyDerivativeByName<Real>(_ns_min_names[i], _eta_name[j]);
162 144 : _d2ns_min[i][j].resize(_neta);
163 :
164 432 : for (unsigned int k = 0; k < _neta; ++k)
165 288 : _d2ns_min[i][j][k] = &getMaterialPropertyDerivativeByName<Real>(
166 288 : _ns_min_names[i], _eta_name[j], _eta_name[k]);
167 : }
168 : }
169 :
170 108 : for (unsigned int i = 0; i < _neta; ++i)
171 : {
172 72 : _eta[i] = &coupledValue("etas", i);
173 72 : _eta_name[i] = coupledName("etas", i);
174 72 : _domegasdeta[i] = &declarePropertyDerivative<Real>("omegas", _eta_name[i]);
175 72 : _d2omegasdetadeta[i].resize(_neta);
176 :
177 180 : for (unsigned int j = 0; j <= i; ++j)
178 108 : _d2omegasdetadeta[j][i] = _d2omegasdetadeta[i][j] =
179 216 : &declarePropertyDerivative<Real>("omegas", _eta_name[j], _eta_name[i]);
180 : }
181 :
182 108 : for (unsigned int i = 0; i < _ndefects; ++i)
183 216 : for (unsigned int j = 0; j < _neta; ++j)
184 144 : _d2omegasdwdeta[i][j] = &declarePropertyDerivative<Real>("omegas", _w_name[i], _eta_name[j]);
185 36 : }
186 :
187 : void
188 441600 : ElectrochemicalSinteringMaterial::computeQpProperties()
189 : {
190 : // Calculate phase switching functions
191 441600 : _hv[_qp] = 0.0;
192 441600 : _dhv[_qp] = 0.0;
193 441600 : _d2hv[_qp] = 0.0;
194 :
195 441600 : if (_phi[_qp] >= 1.0)
196 0 : _hv[_qp] = 1.0;
197 441600 : else if (_phi[_qp] > 0.0)
198 : {
199 0 : _hv[_qp] = _phi[_qp] * _phi[_qp] * _phi[_qp] * (10.0 + _phi[_qp] * (-15.0 + _phi[_qp] * 6.0));
200 0 : _dhv[_qp] = 30.0 * _phi[_qp] * _phi[_qp] * (_phi[_qp] - 1.0) * (_phi[_qp] - 1.0);
201 0 : _d2hv[_qp] = 60.0 * _phi[_qp] * (2.0 * _phi[_qp] - 1.0) * (_phi[_qp] - 1.0);
202 : }
203 :
204 441600 : _hs[_qp] = 1.0 - _hv[_qp];
205 441600 : _dhs[_qp] = -_dhv[_qp];
206 441600 : _d2hs[_qp] = -_d2hv[_qp];
207 :
208 : // Calculate interface switching function
209 441600 : Real phi = _phi[_qp] / _switch;
210 : Real f = 0.0;
211 : Real df = 0.0;
212 : Real d2f = 0.0;
213 441600 : if (phi >= 1.0)
214 : f = 1.0;
215 441600 : else if (phi > 0.0)
216 : {
217 0 : f = phi * phi * phi * (10.0 + phi * (-15.0 + phi * 6.0));
218 0 : df = 30.0 / _switch * phi * phi * (phi - 1.0) * (phi - 1.0);
219 0 : d2f = 60.0 * phi / (_switch * _switch) * (2.0 * phi - 1.0) * (phi - 1.0);
220 : }
221 :
222 : // Permittivity in the void phase (permittivity of free space)
223 441600 : Real eps_0 = 5.52e-2 * _v_scale * _v_scale; // units eV/V^2/nm, change with V_scale if needed
224 :
225 : // Calculate grand potential of void phase
226 : Real sum_omega_v = 0.0;
227 1324800 : for (unsigned int i = 0; i < _ndefects; ++i)
228 883200 : sum_omega_v += -0.5 / (*_kv[i])[_qp] * ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) *
229 883200 : ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) -
230 883200 : _nv_min[i] * ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale);
231 :
232 441600 : _omegav[_qp] = sum_omega_v - 0.5 * eps_0 * _grad_V[_qp] * _grad_V[_qp];
233 :
234 : // Calculate derivatives of grand potential wrt chemical potential
235 1324800 : for (unsigned int i = 0; i < _ndefects; ++i)
236 : {
237 883200 : (*_domegavdw[i])[_qp] =
238 883200 : -((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) / (*_kv[i])[_qp] - _nv_min[i];
239 883200 : (*_d2omegavdw2[i])[_qp] = -1.0 / (*_kv[i])[_qp];
240 : }
241 :
242 : // Calculate solid phase density and potential
243 441600 : switch (_solid_energy)
244 : {
245 : case 0: // PARABOLIC
246 : {
247 : // Calculate grand potential of solid phase and derivatives wrt chemical potentials
248 : Real sum_omega_s = 0.0;
249 0 : for (unsigned int i = 0; i < _ndefects; ++i)
250 : {
251 0 : sum_omega_s += -0.5 / (*_ks[i])[_qp] * ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) *
252 0 : ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) -
253 0 : (*_ns_min[i])[_qp] * ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale);
254 0 : (*_domegasdw[i])[_qp] =
255 0 : -((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) / (*_ks[i])[_qp] -
256 : (*_ns_min[i])[_qp];
257 0 : (*_d2omegasdw2[i])[_qp] = -1.0 / (*_ks[i])[_qp];
258 : }
259 :
260 0 : _omegas[_qp] = sum_omega_s - 0.5 * _eps_r * eps_0 * _grad_V[_qp] * _grad_V[_qp];
261 :
262 : // Calculate derivatives of grand potential wrt order parameters and mixed derivatives
263 0 : for (unsigned int i = 0; i < _neta; ++i)
264 : {
265 : Real sum_domegasdeta = 0.0;
266 0 : for (unsigned int j = 0; j < _ndefects; ++j)
267 : {
268 0 : sum_domegasdeta +=
269 0 : -(*_dns_min[j][i])[_qp] * ((*_w[j])[_qp] - _z[j] * _e * _v[_qp] * _v_scale);
270 0 : (*_d2omegasdwdeta[j][i])[_qp] = -(*_dns_min[j][i])[_qp];
271 : }
272 0 : (*_domegasdeta[i])[_qp] = sum_domegasdeta;
273 :
274 0 : for (unsigned int j = i; j < _neta; ++j)
275 : {
276 : Real sum_d2omegadeta2 = 0.0;
277 0 : for (unsigned int k = 0; k < _ndefects; ++k)
278 0 : sum_d2omegadeta2 +=
279 0 : -(*_d2ns_min[k][i][j])[_qp] * ((*_w[k])[_qp] - _z[k] * _e * _v[_qp] * _v_scale);
280 :
281 0 : (*_d2omegasdetadeta[i][j])[_qp] = (*_d2omegasdetadeta[j][i])[_qp] = sum_d2omegadeta2;
282 : }
283 : }
284 : break;
285 : } // case 0; // PARABOLIC
286 : case 1: // DILUTE
287 : {
288 : // Calculate grand potential of solid phase and derivatives wrt chemical potentials
289 : Real sum_omega_s = 0.0;
290 1324800 : for (unsigned int i = 0; i < _ndefects; ++i)
291 : {
292 883200 : Real n_exp = std::exp(((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
293 883200 : sum_omega_s += -_kB * _temp[_qp] * (*_ns_min[i])[_qp] * n_exp;
294 883200 : (*_domegasdw[i])[_qp] = -(*_ns_min[i])[_qp] * n_exp;
295 883200 : (*_d2omegasdw2[i])[_qp] = -(*_ns_min[i])[_qp] * n_exp / _kB / _temp[_qp];
296 : }
297 :
298 441600 : _omegas[_qp] = sum_omega_s - 0.5 * _eps_r * eps_0 * _grad_V[_qp] * _grad_V[_qp];
299 :
300 : // Calculate derivatives of grand potential wrt order parameters and mixed derivatives
301 1324800 : for (unsigned int i = 0; i < _neta; ++i)
302 : {
303 : Real sum_domegasdeta = 0.0;
304 2649600 : for (unsigned int j = 0; j < _ndefects; ++j)
305 : {
306 : Real n_exp =
307 1766400 : std::exp(((*_w[j])[_qp] - _z[j] * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
308 1766400 : sum_domegasdeta += -_kB / _temp[_qp] * (*_dns_min[j][i])[_qp] * n_exp;
309 1766400 : (*_d2omegasdwdeta[j][i])[_qp] = -(*_dns_min[j][i])[_qp] * n_exp;
310 : }
311 883200 : (*_domegasdeta[i])[_qp] = sum_domegasdeta;
312 :
313 2208000 : for (unsigned int j = i; j < _neta; ++j)
314 : {
315 : Real sum_d2omegadeta2 = 0.0;
316 3974400 : for (unsigned int k = 0; k < _ndefects; ++k)
317 : {
318 : Real n_exp =
319 2649600 : std::exp(((*_w[k])[_qp] - _z[k] * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
320 2649600 : sum_d2omegadeta2 += -_kB / _temp[_qp] * (*_d2ns_min[k][i][j])[_qp] * n_exp;
321 : }
322 :
323 1324800 : (*_d2omegasdetadeta[i][j])[_qp] = (*_d2omegasdetadeta[j][i])[_qp] = sum_d2omegadeta2;
324 : }
325 : }
326 :
327 : break;
328 : } // case 1: // DILUTE
329 0 : case 2: // IDEAL
330 : {
331 0 : mooseError(
332 : "Ideal solution in solid is not yet supported in ElectrochemicalSinteringMaterial");
333 : break;
334 : } // case 2: // IDEAL
335 : } // switch (_solid_energy)
336 :
337 : // thermodynamic parameters
338 441600 : _mu[_qp] = _mu_gb + (_mu_s - _mu_gb) * f;
339 441600 : _kappa[_qp] = _kappa_gb + (_kappa_s - _kappa_gb) * f;
340 441600 : _dmu[_qp] = (_mu_s - _mu_gb) * df;
341 441600 : _dkappa[_qp] = (_kappa_s - _kappa_gb) * df;
342 441600 : _d2mu[_qp] = (_mu_s - _mu_gb) * d2f;
343 441600 : _d2kappa[_qp] = (_kappa_s - _kappa_gb) * d2f;
344 441600 : _gamma[_qp] = 1.5;
345 441600 : } // void ElectrochemicalSinteringMaterial::computeQpProperties()
|
