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 "ElectrochemicalDefectMaterial.h"
11 : #include "libmesh/quadrature.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("PhaseFieldApp", ElectrochemicalDefectMaterial);
15 :
16 : InputParameters
17 94 : ElectrochemicalDefectMaterial::validParams()
18 : {
19 94 : InputParameters params = Material::validParams();
20 94 : params.addClassDescription(
21 : "Calculates density, susceptibility, and derivatives for a defect species in the grand "
22 : "potential sintering model coupled with electrochemistry");
23 188 : params.addRequiredCoupledVarWithAutoBuild(
24 : "etas", "var_name_base", "op_num", "Array of order parameters that describe solid phase");
25 188 : params.addRequiredCoupledVar("chemical_potential", "The name of the chemical potential variable");
26 188 : params.addRequiredCoupledVar("void_op", "The name of the void phase order parameter");
27 188 : params.addRequiredCoupledVar("Temperature", "Name of the temperature variable with units of K");
28 188 : params.addRequiredCoupledVar("electric_potential",
29 : "Name of the electric potential variable with units of V");
30 188 : params.addParam<MaterialPropertyName>(
31 : "void_density_name",
32 : "nv",
33 : "Name of material property to be created for defect number density in the void phase.");
34 188 : params.addParam<MaterialPropertyName>(
35 : "solid_density_name",
36 : "ns",
37 : "Name of material property to be created for defect number density in the solid phase.");
38 188 : params.addParam<MaterialPropertyName>(
39 : "chi_name", "chi", "Name of material property to be created defect susceptibility.");
40 188 : params.addParam<MaterialPropertyName>("solid_energy_coefficient",
41 188 : 1.0,
42 : "Parabolic solid energy coefficient (energy*volume) for "
43 : "defect species. Only used for parabolic energy.");
44 188 : params.addRequiredParam<MaterialPropertyName>(
45 : "void_energy_coefficient",
46 : "Parabolic void energy coefficient (energy*volume) for defect species.");
47 188 : params.addRequiredParam<MaterialPropertyName>(
48 : "min_vacancy_concentration_solid",
49 : "Name of material that determines the minimum in energy wrt defect concentration in the "
50 : "solid phase");
51 188 : params.addRequiredParam<Real>("min_vacancy_concentration_void",
52 : "Minimum in energy wrt defect concentration in the void phase");
53 188 : MooseEnum solid_energy_model("PARABOLIC DILUTE IDEAL", "PARABOLIC");
54 188 : params.addParam<MooseEnum>("solid_energy_model",
55 : solid_energy_model,
56 : "Type of energy function to use for the solid phase.");
57 188 : params.addRequiredParam<int>("defect_charge", "Effective charge of defect species");
58 188 : params.addRequiredParam<Real>("solid_relative_permittivity",
59 : "Solid phase relative permittivity (dimensionless)");
60 188 : params.addParam<Real>("voltage_scale", 1, "Voltage scale (default is for voltage in V)");
61 94 : return params;
62 94 : }
63 :
64 72 : ElectrochemicalDefectMaterial::ElectrochemicalDefectMaterial(const InputParameters & parameters)
65 : : DerivativeMaterialInterface<Material>(parameters),
66 72 : _neta(coupledComponents("etas")),
67 72 : _eta(_neta),
68 72 : _eta_name(_neta),
69 72 : _w(coupledValue("chemical_potential")),
70 72 : _w_name(coupledName("chemical_potential", 0)),
71 72 : _phi(coupledValue("void_op")),
72 72 : _phi_name(coupledName("void_op", 0)),
73 72 : _ns_min_name(getParam<MaterialPropertyName>("min_vacancy_concentration_solid")),
74 72 : _ns_min(getMaterialProperty<Real>(_ns_min_name)),
75 72 : _dns_min(_neta),
76 72 : _d2ns_min(_neta),
77 72 : _temp(coupledValue("Temperature")),
78 72 : _v(coupledValue("electric_potential")),
79 72 : _grad_V(coupledGradient("electric_potential")),
80 144 : _kv(getMaterialProperty<Real>("void_energy_coefficient")),
81 144 : _ks(getMaterialProperty<Real>("solid_energy_coefficient")),
82 144 : _hv(getMaterialPropertyByName<Real>("hv")),
83 144 : _dhv(getMaterialPropertyDerivativeByName<Real>("hv", _phi_name)),
84 144 : _d2hv(getMaterialPropertyDerivativeByName<Real>("hv", _phi_name, _phi_name)),
85 144 : _hs(getMaterialPropertyByName<Real>("hs")),
86 144 : _dhs(getMaterialPropertyDerivativeByName<Real>("hs", _phi_name)),
87 144 : _d2hs(getMaterialPropertyDerivativeByName<Real>("hs", _phi_name, _phi_name)),
88 72 : _chi_name(getParam<MaterialPropertyName>("chi_name")),
89 72 : _chi(declareProperty<Real>(_chi_name)),
90 72 : _dchidphi(declarePropertyDerivative<Real>(_chi_name, _phi_name)),
91 72 : _dchidw(declarePropertyDerivative<Real>(_chi_name, _w_name)),
92 72 : _d2chidphi2(declarePropertyDerivative<Real>(_chi_name, _phi_name, _phi_name)),
93 72 : _d2chidw2(declarePropertyDerivative<Real>(_chi_name, _w_name, _w_name)),
94 72 : _d2chidphidw(declarePropertyDerivative<Real>(_chi_name, _phi_name, _w_name)),
95 72 : _nv_name(getParam<MaterialPropertyName>("void_density_name")),
96 72 : _nv(declareProperty<Real>(_nv_name)),
97 72 : _dnvdw(declarePropertyDerivative<Real>(_nv_name, _w_name)),
98 72 : _ns_name(getParam<MaterialPropertyName>("solid_density_name")),
99 72 : _ns(declareProperty<Real>(_ns_name)),
100 72 : _dnsdw(declarePropertyDerivative<Real>(_ns_name, _w_name)),
101 72 : _d2nsdw2(declarePropertyDerivative<Real>(_ns_name, _w_name, _w_name)),
102 72 : _dns(_neta),
103 72 : _d2nsdwdeta(_neta),
104 72 : _d2ns(_neta),
105 144 : _solid_energy(getParam<MooseEnum>("solid_energy_model")),
106 72 : _kB(8.617343e-5), // eV/K
107 72 : _e(1.0), // to put energy units in eV
108 144 : _z(getParam<int>("defect_charge")),
109 144 : _v_scale(getParam<Real>("voltage_scale")),
110 144 : _eps_r(getParam<Real>("solid_relative_permittivity")),
111 216 : _nv_min(getParam<Real>("min_vacancy_concentration_void"))
112 : {
113 :
114 216 : for (unsigned int i = 0; i < _neta; ++i)
115 : {
116 144 : _eta[i] = &coupledValue("etas", i);
117 144 : _eta_name[i] = coupledName("etas", i);
118 144 : _dns_min[i] = &getMaterialPropertyDerivativeByName<Real>(_ns_min_name, _eta_name[i]);
119 144 : _d2ns_min[i].resize(_neta);
120 144 : _dns[i] = &declarePropertyDerivative<Real>(_ns_name, _eta_name[i]);
121 144 : _d2ns[i].resize(_neta);
122 144 : _d2nsdwdeta[i] = &declarePropertyDerivative<Real>(_ns_name, _w_name, _eta_name[i]);
123 :
124 360 : for (unsigned int j = 0; j <= i; ++j)
125 : {
126 216 : _d2ns_min[j][i] =
127 216 : &getMaterialPropertyDerivativeByName<Real>(_ns_min_name, _eta_name[j], _eta_name[i]);
128 216 : _d2ns[j][i] = &declarePropertyDerivative<Real>(_ns_name, _eta_name[j], _eta_name[i]);
129 : }
130 : }
131 72 : }
132 :
133 : void
134 883200 : ElectrochemicalDefectMaterial::computeQpProperties()
135 : {
136 : // Calculate the void phase density and derivatives
137 883200 : _nv[_qp] = (_w[_qp] - _z * _e * _v[_qp] * _v_scale) / _kv[_qp] + _nv_min;
138 883200 : _dnvdw[_qp] = 1.0 / _kv[_qp]; // susceptibility
139 :
140 : // Calculate solid phase density and derivatives
141 : Real d3nsdw3 = 0.0;
142 883200 : switch (_solid_energy)
143 : {
144 0 : case 0: // PARABOLIC
145 : {
146 0 : _ns[_qp] = (_w[_qp] - _z * _e * _v[_qp] * _v_scale) / _ks[_qp] + _ns_min[_qp];
147 0 : _dnsdw[_qp] = 1.0 / _ks[_qp];
148 0 : _d2nsdw2[_qp] = 0.0;
149 : d3nsdw3 = 0.0;
150 :
151 0 : for (unsigned int i = 0; i < _neta; ++i)
152 : {
153 0 : (*_dns[i])[_qp] = (*_dns_min[i])[_qp];
154 0 : (*_d2nsdwdeta[i])[_qp] = 0.0;
155 0 : for (unsigned int j = i; j < _neta; ++j)
156 : {
157 0 : (*_d2ns[i][j])[_qp] = (*_d2ns_min[i][j])[_qp];
158 : }
159 : }
160 : break;
161 : } // case 0; // PARABOLIC
162 883200 : case 1: // DILUTE
163 : {
164 883200 : Real n_exp = std::exp((_w[_qp] - _z * _e * _v[_qp] * _v_scale) / _kB / _temp[_qp]);
165 883200 : _ns[_qp] = _ns_min[_qp] * n_exp;
166 883200 : _dnsdw[_qp] = _ns[_qp] / _kB / _temp[_qp];
167 883200 : _d2nsdw2[_qp] = _dnsdw[_qp] / _kB / _temp[_qp];
168 883200 : d3nsdw3 = _d2nsdw2[_qp] / _kB / _temp[_qp];
169 :
170 2649600 : for (unsigned int i = 0; i < _neta; ++i)
171 : {
172 1766400 : (*_dns[i])[_qp] = (*_dns_min[i])[_qp] * n_exp;
173 1766400 : (*_d2nsdwdeta[i])[_qp] = (*_dns_min[i])[_qp] * n_exp / _kB / _temp[_qp];
174 4416000 : for (unsigned int j = i; j < _neta; ++j)
175 : {
176 2649600 : (*_d2ns[i][j])[_qp] = (*_d2ns_min[i][j])[_qp] * n_exp;
177 : }
178 : }
179 : break;
180 : } // case 1: // DILUTE
181 0 : case 2: // IDEAL
182 : {
183 0 : mooseError("Ideal solution in solid is not yet supported in ElectrochemicalDefectMaterial");
184 : break;
185 : } // case 2: // IDEAL
186 : } // switch (_solid_energy)
187 :
188 : // Calculate the susceptibilities
189 883200 : _chi[_qp] = _hs[_qp] * _dnsdw[_qp] + _hv[_qp] * _dnvdw[_qp];
190 883200 : _dchidphi[_qp] = _dhs[_qp] * _dnsdw[_qp] + _dhv[_qp] * _dnvdw[_qp];
191 883200 : _dchidw[_qp] = _hs[_qp] * _d2nsdw2[_qp];
192 883200 : _d2chidphi2[_qp] = _d2hs[_qp] * _dnsdw[_qp] + _d2hv[_qp] * _dnvdw[_qp];
193 883200 : _d2chidw2[_qp] = _hs[_qp] * d3nsdw3;
194 883200 : _d2chidphidw[_qp] = _dhs[_qp] * _d2nsdw2[_qp];
195 883200 : } // void ElectrochemicalDefectMaterial::computeQpProperties()
|
