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