Cantera  3.1.0a1
Loading...
Searching...
No Matches
IonFlow.cpp
Go to the documentation of this file.
1//! @file IonFlow.cpp
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
8#include "cantera/oneD/refine.h"
13#include "cantera/base/global.h"
14
15namespace Cantera
16{
17
18IonFlow::IonFlow(ThermoPhase* ph, size_t nsp, size_t points) :
19 StFlow(ph, nsp, points)
20{
21 // make a local copy of species charge
22 for (size_t k = 0; k < m_nsp; k++) {
23 m_speciesCharge.push_back(m_thermo->charge(k));
24 }
25
26 // Find indices for charge of species
27 for (size_t k = 0; k < m_nsp; k++){
28 if (m_speciesCharge[k] != 0){
29 m_kCharge.push_back(k);
30 } else {
31 m_kNeutral.push_back(k);
32 }
33 }
34
35 // Find the index of electron
36 if (m_thermo->speciesIndex("E") != npos ) {
37 m_kElectron = m_thermo->speciesIndex("E");
38 }
39
40 // no bound for electric potential
41 setBounds(c_offset_E, -1.0e20, 1.0e20);
42
43 // Set tighter negative species limit on charged species to avoid
44 // instabilities. Tolerance on electrons is even tighter to account for the
45 // low "molecular" weight.
46 for (size_t k : m_kCharge) {
47 setBounds(c_offset_Y + k, -1e-14, 1.0);
48 }
49 setBounds(c_offset_Y + m_kElectron, -1e-18, 1.0);
50
51 m_refiner->setActive(c_offset_E, false);
52 m_mobility.resize(m_nsp*m_points);
53 m_do_electric_field.resize(m_points,false);
54}
55
56IonFlow::IonFlow(shared_ptr<Solution> sol, const string& id, size_t points)
57 : IonFlow(sol->thermo().get(), sol->thermo()->nSpecies(), points)
58{
59 m_solution = sol;
60 m_id = id;
61 m_kin = m_solution->kinetics().get();
62 m_trans = m_solution->transport().get();
63 if (m_trans->transportModel() == "none") {
64 throw CanteraError("IonFlow::IonFlow",
65 "An appropriate transport model\nshould be set when instantiating the "
66 "Solution ('gas') object.");
67 }
68 m_solution->registerChangedCallback(this, [this]() {
69 setKinetics(m_solution->kinetics());
70 setTransport(m_solution->transport());
71 });
72}
73
74string IonFlow::domainType() const {
75 if (m_isFree) {
76 return "free-ion-flow";
77 }
78 if (m_usesLambda) {
79 return "axisymmetric-ion-flow";
80 }
81 return "unstrained-ion-flow";
82}
83
84void IonFlow::resize(size_t components, size_t points){
85 StFlow::resize(components, points);
87 m_do_species.resize(m_nsp,true);
88 m_do_electric_field.resize(m_points,false);
89}
90
91bool IonFlow::componentActive(size_t n) const
92{
93 if (n == c_offset_E) {
94 return true;
95 } else {
97 }
98}
99
100void IonFlow::updateTransport(double* x, size_t j0, size_t j1)
101{
103 for (size_t j = j0; j < j1; j++) {
104 setGasAtMidpoint(x,j);
105 m_trans->getMobilities(&m_mobility[j*m_nsp]);
107 size_t k = m_kElectron;
108 double tlog = log(m_thermo->temperature());
109 m_mobility[k+m_nsp*j] = poly5(tlog, m_mobi_e_fix.data());
110 double rho = m_thermo->density();
111 double wtm = m_thermo->meanMolecularWeight();
112 m_diff[k+m_nsp*j] = m_wt[k]*rho*poly5(tlog, m_diff_e_fix.data())/wtm;
113 }
114 }
115}
116
117void IonFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
118{
119 if (m_stage == 1) {
120 frozenIonMethod(x,j0,j1);
121 }
122 if (m_stage == 2) {
123 electricFieldMethod(x,j0,j1);
124 }
125}
126
127void IonFlow::frozenIonMethod(const double* x, size_t j0, size_t j1)
128{
129 for (size_t j = j0; j < j1; j++) {
130 double dz = z(j+1) - z(j);
131 double sum = 0.0;
132 for (size_t k : m_kNeutral) {
133 m_flux(k,j) = m_diff[k+m_nsp*j];
134 m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
135 sum -= m_flux(k,j);
136 }
137
138 // correction flux to insure that \sum_k Y_k V_k = 0.
139 for (size_t k : m_kNeutral) {
140 m_flux(k,j) += sum*Y(x,k,j);
141 }
142
143 // flux for ions
144 // Set flux to zero to prevent some fast charged species (such electrons)
145 // to run away
146 for (size_t k : m_kCharge) {
147 m_flux(k,j) = 0;
148 }
149 }
150}
151
152void IonFlow::electricFieldMethod(const double* x, size_t j0, size_t j1)
153{
154 for (size_t j = j0; j < j1; j++) {
155 double rho = density(j);
156 double dz = z(j+1) - z(j);
157
158 // mixture-average diffusion
159 for (size_t k = 0; k < m_nsp; k++) {
160 m_flux(k,j) = m_diff[k+m_nsp*j];
161 m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
162 }
163
164 // ambipolar diffusion
165 double E_ambi = E(x,j);
166 for (size_t k : m_kCharge) {
167 double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1));
168 double drift = rho * Yav * E_ambi
169 * m_speciesCharge[k] * m_mobility[k+m_nsp*j];
170 m_flux(k,j) += drift;
171 }
172
173 // correction flux
174 double sum_flux = 0.0;
175 for (size_t k = 0; k < m_nsp; k++) {
176 sum_flux -= m_flux(k,j); // total net flux
177 }
178 double sum_ion = 0.0;
179 for (size_t k : m_kCharge) {
180 sum_ion += Y(x,k,j);
181 }
182 // The portion of correction for ions is taken off
183 for (size_t k : m_kNeutral) {
184 m_flux(k,j) += Y(x,k,j) / (1-sum_ion) * sum_flux;
185 }
186 }
187}
188
189void IonFlow::setSolvingStage(const size_t stage)
190{
191 if (stage == 1 || stage == 2) {
192 m_stage = stage;
193 } else {
194 throw CanteraError("IonFlow::setSolvingStage",
195 "solution stage must be set to: "
196 "1) frozenIonMethod, "
197 "2) electricFieldEqnMethod");
198 }
199}
200
201//! Evaluate the electric field equation residual
202void IonFlow::evalElectricField(double* x, double* rsd, int* diag,
203 double rdt, size_t jmin, size_t jmax)
204{
205 StFlow::evalElectricField(x, rsd, diag, rdt, jmin, jmax);
206 if (m_stage != 2) {
207 return;
208 }
209
210 if (jmin == 0) { // left boundary
211 rsd[index(c_offset_E, jmin)] = E(x,jmin);
212 }
213
214 if (jmax == m_points - 1) { // right boundary
215 rsd[index(c_offset_E, jmax)] = dEdz(x,jmax) - rho_e(x,jmax) / epsilon_0;
216 }
217
218 // j0 and j1 are constrained to only interior points
219 size_t j0 = std::max<size_t>(jmin, 1);
220 size_t j1 = std::min(jmax, m_points - 2);
221 for (size_t j = j0; j <= j1; j++) {
222 rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
223 diag[index(c_offset_E, j)] = 0;
224 }
225}
226
227void IonFlow::evalSpecies(double* x, double* rsd, int* diag,
228 double rdt, size_t jmin, size_t jmax)
229{
230 StFlow::evalSpecies(x, rsd, diag, rdt, jmin, jmax);
231 if (m_stage != 2) {
232 return;
233 }
234
235 if (jmin == 0) { // left boundary
236 // enforcing the flux for charged species is difficult
237 // since charged species are also affected by electric
238 // force, so Neumann boundary condition is used.
239 for (size_t k : m_kCharge) {
240 rsd[index(c_offset_Y + k, jmin)] = Y(x,k,jmin) - Y(x,k,jmin + 1);
241 }
242 }
243}
244
246{
247 bool changed = false;
248 if (j == npos) {
249 for (size_t i = 0; i < m_points; i++) {
250 if (!m_do_electric_field[i]) {
251 changed = true;
252 }
253 m_do_electric_field[i] = true;
254 }
255 } else {
256 if (!m_do_electric_field[j]) {
257 changed = true;
258 }
259 m_do_electric_field[j] = true;
260 }
261 m_refiner->setActive(c_offset_U, true);
262 m_refiner->setActive(c_offset_V, true);
263 m_refiner->setActive(c_offset_T, true);
264 m_refiner->setActive(c_offset_E, true);
265 if (changed) {
267 }
268}
269
271{
272 bool changed = false;
273 if (j == npos) {
274 for (size_t i = 0; i < m_points; i++) {
275 if (m_do_electric_field[i]) {
276 changed = true;
277 }
278 m_do_electric_field[i] = false;
279 }
280 } else {
281 if (m_do_electric_field[j]) {
282 changed = true;
283 }
284 m_do_electric_field[j] = false;
285 }
286 m_refiner->setActive(c_offset_U, false);
287 m_refiner->setActive(c_offset_V, false);
288 m_refiner->setActive(c_offset_T, false);
289 m_refiner->setActive(c_offset_E, false);
290 if (changed) {
292 }
293}
294
295void IonFlow::setElectronTransport(vector<double>& tfix, vector<double>& diff_e,
296 vector<double>& mobi_e)
297{
299 size_t degree = 5;
300 size_t n = tfix.size();
301 vector<double> tlog;
302 for (size_t i = 0; i < n; i++) {
303 tlog.push_back(log(tfix[i]));
304 }
305 vector<double> w(n, -1.0);
306 m_diff_e_fix.resize(degree + 1);
307 m_mobi_e_fix.resize(degree + 1);
308 polyfit(n, degree, tlog.data(), diff_e.data(), w.data(), m_diff_e_fix.data());
309 polyfit(n, degree, tlog.data(), mobi_e.data(), w.data(), m_mobi_e_fix.data());
310}
311
312void IonFlow::_finalize(const double* x)
313{
315
316 bool p = m_do_electric_field[0];
317 if (p) {
319 }
320}
321
322}
Headers for the Transport object, which is the virtual base class for all transport property evaluato...
Base class for exceptions thrown by Cantera classes.
shared_ptr< Solution > m_solution
Composite thermo/kinetics/transport handler.
Definition Domain1D.h:567
size_t m_points
Number of grid points.
Definition Domain1D.h:537
string m_id
Identity tag for the domain.
Definition Domain1D.h:560
void needJacUpdate()
Set this if something has changed in the governing equations (for example, the value of a constant ha...
Definition Domain1D.cpp:95
This class models the ion transportation in a flame.
Definition IonFlow.h:29
vector< size_t > m_kCharge
index of species with charges
Definition IonFlow.h:125
void electricFieldMethod(const double *x, size_t j0, size_t j1)
Solving phase two: the electric field equation is added coupled by the electrical drift.
Definition IonFlow.cpp:152
double E(const double *x, size_t j) const
electric field
Definition IonFlow.h:144
vector< bool > m_do_electric_field
flag for solving electric field or not
Definition IonFlow.h:116
size_t m_kElectron
index of electron
Definition IonFlow.h:141
void frozenIonMethod(const double *x, size_t j0, size_t j1)
Solving phase one: the fluxes of charged species are turned off.
Definition IonFlow.cpp:127
void resize(size_t components, size_t points) override
Resize the domain to have nv components and np grid points.
Definition IonFlow.cpp:84
void setElectronTransport(vector< double > &tfix, vector< double > &diff_e, vector< double > &mobi_e)
Sometimes it is desired to carry out the simulation using a specified electron transport profile,...
Definition IonFlow.cpp:295
void evalElectricField(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax) override
Evaluate the electric field equation residual by Gauss's law.
Definition IonFlow.cpp:202
double rho_e(double *x, size_t j) const
total charge density
Definition IonFlow.h:158
void updateTransport(double *x, size_t j0, size_t j1) override
Update the transport properties at grid points in the range from j0 to j1, based on solution x.
Definition IonFlow.cpp:100
vector< double > m_mobility
mobility
Definition IonFlow.h:135
void updateDiffFluxes(const double *x, size_t j0, size_t j1) override
Update the diffusive mass fluxes.
Definition IonFlow.cpp:117
void evalSpecies(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax) override
Evaluate the species equations' residual.
Definition IonFlow.cpp:227
size_t m_stage
solving stage
Definition IonFlow.h:138
void _finalize(const double *x) override
In some cases, a domain may need to set parameters that depend on the initial solution estimate.
Definition IonFlow.cpp:312
void setSolvingStage(const size_t stage) override
Solving stage mode for handling ionized species (used by IonFlow specialization)
Definition IonFlow.cpp:189
bool m_import_electron_transport
flag for importing transport of electron
Definition IonFlow.h:119
void solveElectricField(size_t j=npos) override
Set to solve electric field in a point (used by IonFlow specialization)
Definition IonFlow.cpp:245
void fixElectricField(size_t j=npos) override
Set to fix voltage in a point (used by IonFlow specialization)
Definition IonFlow.cpp:270
string domainType() const override
Domain type flag.
Definition IonFlow.cpp:74
vector< size_t > m_kNeutral
index of neutral species
Definition IonFlow.h:128
vector< double > m_mobi_e_fix
coefficients of polynomial fitting of fixed electron transport profile
Definition IonFlow.h:131
bool componentActive(size_t n) const override
Returns true if the specified component is an active part of the solver state.
Definition IonFlow.cpp:91
vector< double > m_speciesCharge
electrical properties
Definition IonFlow.h:122
double temperature() const
Temperature (K).
Definition Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:655
virtual double density() const
Density (kg/m^3).
Definition Phase.h:587
void setTransport(shared_ptr< Transport > trans) override
Set transport model to existing instance.
Definition StFlow.cpp:140
void setKinetics(shared_ptr< Kinetics > kin) override
Set the kinetics manager.
Definition StFlow.cpp:134
void resize(size_t components, size_t points) override
Change the grid size. Called after grid refinement.
Definition StFlow.cpp:160
vector< double > m_diff
Array of size m_nsp by m_points for saving density times diffusion coefficient times species molar ma...
Definition StFlow.h:611
virtual bool componentActive(size_t n) const
Returns true if the specified component is an active part of the solver state.
Definition StFlow.cpp:745
virtual void evalSpecies(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the species equations' residuals.
Definition StFlow.cpp:560
virtual void evalElectricField(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the electric field equation residual to be zero everywhere.
Definition StFlow.cpp:599
void _finalize(const double *x) override
In some cases, a domain may need to set parameters that depend on the initial solution estimate.
Definition StFlow.cpp:249
size_t m_nsp
Number of species in the mechanism.
Definition StFlow.h:625
void setGasAtMidpoint(const double *x, size_t j)
Set the gas state to be consistent with the solution at the midpoint between j and j + 1.
Definition StFlow.cpp:237
virtual void updateTransport(double *x, size_t j0, size_t j1)
Update the transport properties at grid points in the range from j0 to j1, based on solution x.
Definition StFlow.cpp:608
virtual string transportModel() const
Identifies the model represented by this Transport object.
Definition Transport.h:93
virtual void getMobilities(double *const mobil_e)
Get the Electrical mobilities (m^2/V/s).
Definition Transport.h:182
Header for a file containing miscellaneous numerical functions.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
Definition utilities.h:141
double polyfit(size_t n, size_t deg, const double *xp, const double *yp, const double *wp, double *pp)
Fits a polynomial function to a set of data points.
Definition polyfit.cpp:14
const double epsilon_0
Permittivity of free space [F/m].
Definition ct_defs.h:137
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
@ c_offset_U
axial velocity
Definition StFlow.h:25
@ c_offset_V
strain rate
Definition StFlow.h:26
@ c_offset_E
electric field equation
Definition StFlow.h:29
@ c_offset_Y
mass fractions
Definition StFlow.h:30
@ c_offset_T
temperature
Definition StFlow.h:27
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...