Cantera  2.4.0
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 http://www.cantera.org/license.txt for license and copyright information.
5 
6 #include "cantera/oneD/IonFlow.h"
7 #include "cantera/oneD/StFlow.h"
8 #include "cantera/base/ctml.h"
10 #include "cantera/numerics/funcs.h"
12 
13 using namespace std;
14 
15 namespace Cantera
16 {
17 
18 IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
19  StFlow(ph, nsp, points),
20  m_import_electron_transport(false),
21  m_stage(1),
22  m_kElectron(npos)
23 {
24  // make a local copy of species charge
25  for (size_t k = 0; k < m_nsp; k++) {
26  m_speciesCharge.push_back(m_thermo->charge(k));
27  }
28 
29  // Find indices for charge of species
30  for (size_t k = 0; k < m_nsp; k++){
31  if (m_speciesCharge[k] != 0){
32  m_kCharge.push_back(k);
33  } else {
34  m_kNeutral.push_back(k);
35  }
36  }
37 
38  // Find the index of electron
39  if (m_thermo->speciesIndex("E") != npos ) {
40  m_kElectron = m_thermo->speciesIndex("E");
41  }
42 
43  // no bound for electric potential
44  setBounds(c_offset_E, -1.0e20, 1.0e20);
45 
46  // Set tighter negative species limit on charged species to avoid
47  // instabilities. Tolerance on electrons is even tighter to account for the
48  // low "molecular" weight.
49  for (size_t k : m_kCharge) {
50  setBounds(c_offset_Y + k, -1e-14, 1.0);
51  }
52  setBounds(c_offset_Y + m_kElectron, -1e-18, 1.0);
53 
54  m_refiner->setActive(c_offset_E, false);
55  m_mobility.resize(m_nsp*m_points);
56  m_do_electric_field.resize(m_points,false);
57 }
58 
59 void IonFlow::resize(size_t components, size_t points){
60  StFlow::resize(components, points);
61  m_mobility.resize(m_nsp*m_points);
62  m_do_species.resize(m_nsp,true);
63  m_do_electric_field.resize(m_points,false);
64 }
65 
66 void IonFlow::updateTransport(double* x, size_t j0, size_t j1)
67 {
68  StFlow::updateTransport(x,j0,j1);
69  for (size_t j = j0; j < j1; j++) {
70  setGasAtMidpoint(x,j);
71  m_trans->getMobilities(&m_mobility[j*m_nsp]);
73  size_t k = m_kElectron;
74  double tlog = log(m_thermo->temperature());
75  m_mobility[k+m_nsp*j] = poly5(tlog, m_mobi_e_fix.data());
76  m_diff[k+m_nsp*j] = poly5(tlog, m_diff_e_fix.data());
77  }
78  }
79 }
80 
81 void IonFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
82 {
83  if (m_stage == 1) {
84  frozenIonMethod(x,j0,j1);
85  }
86  if (m_stage == 2) {
87  electricFieldMethod(x,j0,j1);
88  }
89 }
90 
91 void IonFlow::frozenIonMethod(const double* x, size_t j0, size_t j1)
92 {
93  for (size_t j = j0; j < j1; j++) {
94  double wtm = m_wtm[j];
95  double rho = density(j);
96  double dz = z(j+1) - z(j);
97  double sum = 0.0;
98  for (size_t k : m_kNeutral) {
99  m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
100  m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
101  sum -= m_flux(k,j);
102  }
103 
104  // correction flux to insure that \sum_k Y_k V_k = 0.
105  for (size_t k : m_kNeutral) {
106  m_flux(k,j) += sum*Y(x,k,j);
107  }
108 
109  // flux for ions
110  // Set flux to zero to prevent some fast charged species (e.g. electron)
111  // to run away
112  for (size_t k : m_kCharge) {
113  m_flux(k,j) = 0;
114  }
115  }
116 }
117 
118 void IonFlow::electricFieldMethod(const double* x, size_t j0, size_t j1)
119 {
120  for (size_t j = j0; j < j1; j++) {
121  double wtm = m_wtm[j];
122  double rho = density(j);
123  double dz = z(j+1) - z(j);
124 
125  // mixture-average diffusion
126  double sum = 0.0;
127  for (size_t k = 0; k < m_nsp; k++) {
128  m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm);
129  m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz;
130  sum -= m_flux(k,j);
131  }
132 
133  // ambipolar diffusion
134  double E_ambi = E(x,j);
135  for (size_t k : m_kCharge) {
136  double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1));
137  double drift = rho * Yav * E_ambi
138  * m_speciesCharge[k] * m_mobility[k+m_nsp*j];
139  m_flux(k,j) += drift;
140  }
141 
142  // correction flux
143  double sum_flux = 0.0;
144  for (size_t k = 0; k < m_nsp; k++) {
145  sum_flux -= m_flux(k,j); // total net flux
146  }
147  double sum_ion = 0.0;
148  for (size_t k : m_kCharge) {
149  sum_ion += Y(x,k,j);
150  }
151  // The portion of correction for ions is taken off
152  for (size_t k : m_kNeutral) {
153  m_flux(k,j) += Y(x,k,j) / (1-sum_ion) * sum_flux;
154  }
155  }
156 }
157 
158 void IonFlow::setSolvingStage(const size_t stage)
159 {
160  if (stage == 1 || stage == 2) {
161  m_stage = stage;
162  } else {
163  throw CanteraError("IonFlow::updateDiffFluxes",
164  "solution stage must be set to: "
165  "1) frozenIonMethod, "
166  "2) electricFieldEqnMethod");
167  }
168 }
169 
170 void IonFlow::evalResidual(double* x, double* rsd, int* diag,
171  double rdt, size_t jmin, size_t jmax)
172 {
173  StFlow::evalResidual(x, rsd, diag, rdt, jmin, jmax);
174  if (m_stage != 2) {
175  return;
176  }
177 
178  for (size_t j = jmin; j <= jmax; j++) {
179  if (j == 0) {
180  // enforcing the flux for charged species is difficult
181  // since charged species are also affected by electric
182  // force, so Neumann boundary condition is used.
183  for (size_t k : m_kCharge) {
184  rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1);
185  }
186  rsd[index(c_offset_E, j)] = E(x,0);
187  diag[index(c_offset_E, j)] = 0;
188  } else if (j == m_points - 1) {
189  rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
190  diag[index(c_offset_E, j)] = 0;
191  } else {
192  //-----------------------------------------------
193  // Electric field by Gauss's law
194  //
195  // dE/dz = e/eps_0 * sum(q_k*n_k)
196  //
197  // E = -dV/dz
198  //-----------------------------------------------
199  rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
200  diag[index(c_offset_E, j)] = 0;
201  }
202  }
203 }
204 
206 {
207  bool changed = false;
208  if (j == npos) {
209  for (size_t i = 0; i < m_points; i++) {
210  if (!m_do_electric_field[i]) {
211  changed = true;
212  }
213  m_do_electric_field[i] = true;
214  }
215  } else {
216  if (!m_do_electric_field[j]) {
217  changed = true;
218  }
219  m_do_electric_field[j] = true;
220  }
221  m_refiner->setActive(c_offset_U, true);
222  m_refiner->setActive(c_offset_V, true);
223  m_refiner->setActive(c_offset_T, true);
224  m_refiner->setActive(c_offset_E, true);
225  if (changed) {
226  needJacUpdate();
227  }
228 }
229 
231 {
232  bool changed = false;
233  if (j == npos) {
234  for (size_t i = 0; i < m_points; i++) {
235  if (m_do_electric_field[i]) {
236  changed = true;
237  }
238  m_do_electric_field[i] = false;
239  }
240  } else {
241  if (m_do_electric_field[j]) {
242  changed = true;
243  }
244  m_do_electric_field[j] = false;
245  }
246  m_refiner->setActive(c_offset_U, false);
247  m_refiner->setActive(c_offset_V, false);
248  m_refiner->setActive(c_offset_T, false);
249  m_refiner->setActive(c_offset_E, false);
250  if (changed) {
251  needJacUpdate();
252  }
253 }
254 
256  vector_fp& mobi_e)
257 {
259  size_t degree = 5;
260  size_t n = tfix.size();
261  vector_fp tlog;
262  for (size_t i = 0; i < n; i++) {
263  tlog.push_back(log(tfix[i]));
264  }
265  vector_fp w(n, -1.0);
266  m_diff_e_fix.resize(degree + 1);
267  m_mobi_e_fix.resize(degree + 1);
268  polyfit(n, degree, tlog.data(), diff_e.data(), w.data(), m_diff_e_fix.data());
269  polyfit(n, degree, tlog.data(), mobi_e.data(), w.data(), m_mobi_e_fix.data());
270 }
271 
272 void IonFlow::_finalize(const double* x)
273 {
275 
276  bool p = m_do_electric_field[0];
277  if (p) {
279  }
280 }
281 
282 }
void solveElectricField(size_t j=npos)
set to solve electric field on a point
Definition: IonFlow.cpp:205
virtual void resize(size_t components, size_t points)
Change the grid size. Called after grid refinement.
Definition: IonFlow.cpp:59
vector_fp m_mobility
mobility
Definition: IonFlow.h:98
std::vector< thermo_t * > m_thermo
m_thermo is a vector of pointers to ThermoPhase objects that are involved with this kinetics operator...
Definition: Kinetics.h:882
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: IonFlow.cpp:66
bool m_import_electron_transport
flag for importing transport of electron
Definition: IonFlow.h:82
size_t m_kElectron
index of electron
Definition: IonFlow.h:104
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
Definition: utilities.h:455
virtual void _finalize(const doublereal *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate...
Definition: StFlow.cpp:186
Headers for the Transport object, which is the virtual base class for all transport property evaluato...
STL namespace.
double E(const double *x, size_t j) const
electric field
Definition: IonFlow.h:107
std::vector< bool > m_do_electric_field
flag for solving electric field or not
Definition: IonFlow.h:79
virtual void evalResidual(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Evaluate the residual function.
Definition: StFlow.cpp:290
double rho_e(double *x, size_t j) const
total charge density
Definition: IonFlow.h:121
std::vector< size_t > m_kNeutral
index of neutral species
Definition: IonFlow.h:91
virtual void evalResidual(double *x, double *rsd, int *diag, double rdt, size_t jmin, size_t jmax)
Definition: IonFlow.cpp:170
vector_int m_speciesCharge
electrical properties
Definition: IonFlow.h:85
virtual void _finalize(const double *x)
In some cases, a domain may need to set parameters that depend on the initial solution estimate...
Definition: IonFlow.cpp:272
virtual void updateDiffFluxes(const double *x, size_t j0, size_t j1)
Update the diffusive mass fluxes.
Definition: IonFlow.cpp:81
virtual void setSolvingStage(const size_t phase)
set the solving stage
Definition: IonFlow.cpp:158
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void getMobilities(doublereal *const mobil_e)
Get the Electrical mobilities (m^2/V/s).
virtual 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:118
int m_stage
solving stage
Definition: IonFlow.h:101
virtual void updateTransport(doublereal *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:484
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
void setGasAtMidpoint(const doublereal *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:174
std::vector< size_t > m_kCharge
index of species with charges
Definition: IonFlow.h:88
const doublereal epsilon_0
Permittivity of free space in F/m.
Definition: ct_defs.h:106
void fixElectricField(size_t j=npos)
set to fix voltage on a point
Definition: IonFlow.cpp:230
vector_fp m_mobi_e_fix
coefficients of polynomial fitting of fixed electron transport profile
Definition: IonFlow.h:94
void needJacUpdate()
Definition: Domain1D.cpp:102
virtual void resize(size_t components, size_t points)
Change the grid size. Called after grid refinement.
Definition: StFlow.cpp:97
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
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
Header for a file containing miscellaneous numerical functions.
void setElectronTransport(vector_fp &tfix, vector_fp &diff_e, vector_fp &mobi_e)
Sometimes it is desired to carry out the simulation using a specified electron transport profile...
Definition: IonFlow.cpp:255
virtual 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:91