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 31 : ElectrochemicalSinteringMaterial::validParams()
18 : {
19 31 : InputParameters params = Material::validParams();
20 31 : params.addClassDescription("Includes switching and thermodynamic properties for the grand "
21 : "potential sintering model coupled with electrochemistry");
22 62 : params.addRequiredCoupledVarWithAutoBuild(
23 : "etas", "var_name_base", "op_num", "Array of order parameters that describe solid phase");
24 62 : params.addRequiredCoupledVar("chemical_potentials",
25 : "The name of the chemical potential variables for defects");
26 62 : params.addRequiredCoupledVar("void_op", "The name of the void phase order parameter");
27 62 : params.addRequiredCoupledVar("Temperature", "Name of the temperature variable with units of K");
28 62 : params.addRequiredCoupledVar("electric_potential",
29 : "Name of the electric potential variable with units of V");
30 62 : 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 62 : 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 62 : params.addParam<Real>("surface_energy", 19.7, "Surface energy in units of problem (energy/area)");
40 62 : params.addParam<Real>(
41 62 : "grainboundary_energy", 9.86, "Grain boundary energy in units of problem (energy/area)");
42 62 : params.addParam<Real>("int_width", 1, "Interface width in units of problem (length)");
43 93 : params.addRangeCheckedParam<Real>(
44 : "surface_switch_value",
45 62 : 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 62 : 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 62 : 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 62 : MooseEnum solid_energy_model("PARABOLIC DILUTE IDEAL", "PARABOLIC");
59 62 : params.addParam<MooseEnum>("solid_energy_model",
60 : solid_energy_model,
61 : "Type of energy function to use for the solid phase.");
62 62 : params.addRequiredParam<std::vector<int>>(
63 : "defect_charges",
64 : "Vector of charges of defect species. Place in same order as chemical_potentials!");
65 62 : params.addRequiredParam<Real>("solid_relative_permittivity",
66 : "Solid phase relative permittivity (dimensionless)");
67 62 : params.addParam<Real>("voltage_scale", 1, "Voltage scale (default is for voltage in V)");
68 31 : return params;
69 31 : }
70 :
71 24 : ElectrochemicalSinteringMaterial::ElectrochemicalSinteringMaterial(
72 24 : const InputParameters & parameters)
73 : : DerivativeMaterialInterface<Material>(parameters),
74 24 : _neta(coupledComponents("etas")),
75 24 : _eta(_neta),
76 24 : _eta_name(_neta),
77 24 : _ndefects(coupledComponents("chemical_potentials")),
78 24 : _w(coupledValues("chemical_potentials")),
79 24 : _w_name(_ndefects),
80 24 : _phi(coupledValue("void_op")),
81 24 : _phi_name(coupledName("void_op", 0)),
82 48 : _ns_min_names(getParam<std::vector<MaterialPropertyName>>("min_vacancy_concentrations_solid")),
83 24 : _ns_min(_ndefects),
84 24 : _dns_min(_ndefects),
85 24 : _d2ns_min(_ndefects),
86 24 : _temp(coupledValue("Temperature")),
87 24 : _v(coupledValue("electric_potential")),
88 24 : _grad_V(coupledGradient("electric_potential")),
89 48 : _kv_names(getParam<std::vector<MaterialPropertyName>>("void_energy_coefficients")),
90 24 : _kv(_ndefects),
91 48 : _ks_names(getParam<std::vector<MaterialPropertyName>>("solid_energy_coefficients")),
92 24 : _ks(_ndefects),
93 24 : _hv(declareProperty<Real>("hv")),
94 24 : _dhv(declarePropertyDerivative<Real>("hv", _phi_name)),
95 24 : _d2hv(declarePropertyDerivative<Real>("hv", _phi_name, _phi_name)),
96 24 : _hs(declareProperty<Real>("hs")),
97 24 : _dhs(declarePropertyDerivative<Real>("hs", _phi_name)),
98 24 : _d2hs(declarePropertyDerivative<Real>("hs", _phi_name, _phi_name)),
99 24 : _omegav(declareProperty<Real>("omegav")),
100 24 : _domegavdw(_ndefects),
101 24 : _d2omegavdw2(_ndefects),
102 24 : _omegas(declareProperty<Real>("omegas")),
103 24 : _domegasdw(_ndefects),
104 24 : _d2omegasdw2(_ndefects),
105 24 : _domegasdeta(_neta),
106 24 : _d2omegasdwdeta(_ndefects),
107 24 : _d2omegasdetadeta(_neta),
108 24 : _mu(declareProperty<Real>("mu")),
109 24 : _dmu(declarePropertyDerivative<Real>("mu", _phi_name)),
110 24 : _d2mu(declarePropertyDerivative<Real>("mu", _phi_name, _phi_name)),
111 24 : _kappa(declareProperty<Real>("kappa")),
112 24 : _dkappa(declarePropertyDerivative<Real>("kappa", _phi_name)),
113 24 : _d2kappa(declarePropertyDerivative<Real>("kappa", _phi_name, _phi_name)),
114 24 : _gamma(declareProperty<Real>("gamma")),
115 48 : _sigma_s(getParam<Real>("surface_energy")),
116 48 : _sigma_gb(getParam<Real>("grainboundary_energy")),
117 48 : _int_width(getParam<Real>("int_width")),
118 48 : _switch(getParam<Real>("surface_switch_value")),
119 48 : _solid_energy(getParam<MooseEnum>("solid_energy_model")),
120 24 : _mu_s(6.0 * _sigma_s / _int_width),
121 24 : _mu_gb(6.0 * _sigma_gb / _int_width),
122 24 : _kappa_s(0.75 * _sigma_s * _int_width),
123 24 : _kappa_gb(0.75 * _sigma_gb * _int_width),
124 24 : _kB(8.617343e-5), // eV/K
125 24 : _e(1.0), // to put energy units in eV
126 48 : _z(getParam<std::vector<int>>("defect_charges")),
127 48 : _v_scale(getParam<Real>("voltage_scale")),
128 48 : _eps_r(getParam<Real>("solid_relative_permittivity")),
129 72 : _nv_min(getParam<std::vector<Real>>("min_vacancy_concentrations_void"))
130 : {
131 24 : if (_ndefects != _z.size())
132 0 : paramError("defect_charges", "Need to pass in as many defect_charges as chemical_potentials");
133 24 : 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 24 : 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 24 : 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 24 : 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 72 : for (unsigned int i = 0; i < _ndefects; ++i)
147 : {
148 48 : _w_name[i] = coupledName("chemical_potentials", i);
149 48 : _ns_min[i] = &getMaterialPropertyByName<Real>(_ns_min_names[i]);
150 48 : _dns_min[i].resize(_neta);
151 48 : _d2ns_min[i].resize(_neta);
152 48 : _d2omegasdwdeta[i].resize(_neta);
153 48 : _kv[i] = &getMaterialPropertyByName<Real>(_kv_names[i]);
154 48 : _ks[i] = &getMaterialPropertyByName<Real>(_ks_names[i]);
155 48 : _domegavdw[i] = &declarePropertyDerivative<Real>("omegav", _w_name[i]);
156 48 : _d2omegavdw2[i] = &declarePropertyDerivative<Real>("omegav", _w_name[i], _w_name[i]);
157 48 : _domegasdw[i] = &declarePropertyDerivative<Real>("omegas", _w_name[i]);
158 48 : _d2omegasdw2[i] = &declarePropertyDerivative<Real>("omegas", _w_name[i], _w_name[i]);
159 144 : for (unsigned int j = 0; j < _neta; ++j)
160 : {
161 96 : _dns_min[i][j] = &getMaterialPropertyDerivativeByName<Real>(_ns_min_names[i], _eta_name[j]);
162 96 : _d2ns_min[i][j].resize(_neta);
163 :
164 288 : for (unsigned int k = 0; k < _neta; ++k)
165 192 : _d2ns_min[i][j][k] = &getMaterialPropertyDerivativeByName<Real>(
166 192 : _ns_min_names[i], _eta_name[j], _eta_name[k]);
167 : }
168 : }
169 :
170 72 : for (unsigned int i = 0; i < _neta; ++i)
171 : {
172 48 : _eta[i] = &coupledValue("etas", i);
173 48 : _eta_name[i] = coupledName("etas", i);
174 48 : _domegasdeta[i] = &declarePropertyDerivative<Real>("omegas", _eta_name[i]);
175 48 : _d2omegasdetadeta[i].resize(_neta);
176 :
177 120 : for (unsigned int j = 0; j <= i; ++j)
178 72 : _d2omegasdetadeta[j][i] = _d2omegasdetadeta[i][j] =
179 144 : &declarePropertyDerivative<Real>("omegas", _eta_name[j], _eta_name[i]);
180 : }
181 :
182 72 : for (unsigned int i = 0; i < _ndefects; ++i)
183 144 : for (unsigned int j = 0; j < _neta; ++j)
184 96 : _d2omegasdwdeta[i][j] = &declarePropertyDerivative<Real>("omegas", _w_name[i], _eta_name[j]);
185 24 : }
186 :
187 : void
188 353600 : ElectrochemicalSinteringMaterial::computeQpProperties()
189 : {
190 : // Calculate phase switching functions
191 353600 : _hv[_qp] = 0.0;
192 353600 : _dhv[_qp] = 0.0;
193 353600 : _d2hv[_qp] = 0.0;
194 :
195 353600 : if (_phi[_qp] >= 1.0)
196 0 : _hv[_qp] = 1.0;
197 353600 : 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 353600 : _hs[_qp] = 1.0 - _hv[_qp];
205 353600 : _dhs[_qp] = -_dhv[_qp];
206 353600 : _d2hs[_qp] = -_d2hv[_qp];
207 :
208 : // Calculate interface switching function
209 353600 : Real phi = _phi[_qp] / _switch;
210 : Real f = 0.0;
211 : Real df = 0.0;
212 : Real d2f = 0.0;
213 353600 : if (phi >= 1.0)
214 : f = 1.0;
215 353600 : 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 353600 : 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 1060800 : for (unsigned int i = 0; i < _ndefects; ++i)
228 707200 : sum_omega_v += -0.5 / (*_kv[i])[_qp] * ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) *
229 707200 : ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) -
230 707200 : _nv_min[i] * ((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale);
231 :
232 353600 : _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 1060800 : for (unsigned int i = 0; i < _ndefects; ++i)
236 : {
237 707200 : (*_domegavdw[i])[_qp] =
238 707200 : -((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) / (*_kv[i])[_qp] - _nv_min[i];
239 707200 : (*_d2omegavdw2[i])[_qp] = -1.0 / (*_kv[i])[_qp];
240 : }
241 :
242 : // Calculate solid phase density and potential
243 353600 : 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 1060800 : for (unsigned int i = 0; i < _ndefects; ++i)
291 : {
292 707200 : Real n_exp = std::exp(((*_w[i])[_qp] - _z[i] * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
293 707200 : sum_omega_s += -_kB * _temp[_qp] * (*_ns_min[i])[_qp] * n_exp;
294 707200 : (*_domegasdw[i])[_qp] = -(*_ns_min[i])[_qp] * n_exp;
295 707200 : (*_d2omegasdw2[i])[_qp] = -(*_ns_min[i])[_qp] * n_exp / _kB / _temp[_qp];
296 : }
297 :
298 353600 : _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 1060800 : for (unsigned int i = 0; i < _neta; ++i)
302 : {
303 : Real sum_domegasdeta = 0.0;
304 2121600 : for (unsigned int j = 0; j < _ndefects; ++j)
305 : {
306 : Real n_exp =
307 1414400 : std::exp(((*_w[j])[_qp] - _z[j] * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
308 1414400 : sum_domegasdeta += -_kB / _temp[_qp] * (*_dns_min[j][i])[_qp] * n_exp;
309 1414400 : (*_d2omegasdwdeta[j][i])[_qp] = -(*_dns_min[j][i])[_qp] * n_exp;
310 : }
311 707200 : (*_domegasdeta[i])[_qp] = sum_domegasdeta;
312 :
313 1768000 : for (unsigned int j = i; j < _neta; ++j)
314 : {
315 : Real sum_d2omegadeta2 = 0.0;
316 3182400 : for (unsigned int k = 0; k < _ndefects; ++k)
317 : {
318 : Real n_exp =
319 2121600 : std::exp(((*_w[k])[_qp] - _z[k] * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
320 2121600 : sum_d2omegadeta2 += -_kB / _temp[_qp] * (*_d2ns_min[k][i][j])[_qp] * n_exp;
321 : }
322 :
323 1060800 : (*_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 353600 : _mu[_qp] = _mu_gb + (_mu_s - _mu_gb) * f;
339 353600 : _kappa[_qp] = _kappa_gb + (_kappa_s - _kappa_gb) * f;
340 353600 : _dmu[_qp] = (_mu_s - _mu_gb) * df;
341 353600 : _dkappa[_qp] = (_kappa_s - _kappa_gb) * df;
342 353600 : _d2mu[_qp] = (_mu_s - _mu_gb) * d2f;
343 353600 : _d2kappa[_qp] = (_kappa_s - _kappa_gb) * d2f;
344 353600 : _gamma[_qp] = 1.5;
345 353600 : } // void ElectrochemicalSinteringMaterial::computeQpProperties()
|
