Memosa-FVM  0.2
Public Types | Public Member Functions | Private Attributes | List of all members
LinearizeSpeciesInterface< X, Diag, OffDiag > Class Template Reference

#include <LinearizeSpeciesInterface.h>

Collaboration diagram for LinearizeSpeciesInterface< X, Diag, OffDiag >:
Collaboration graph
[legend]

Public Types

typedef NumTypeTraits< X >
::T_Scalar 		T_Scalar
 
typedef CRMatrix< Diag,
OffDiag, X > 		CCMatrix
 
typedef CCMatrix::DiagArray DiagArray
 
typedef Array< T_Scalar > TArray
 
typedef Array< X > XArray
 
typedef Vector< T_Scalar, 3 > VectorT3
 
typedef Array< VectorT3 > VectorT3Array
 

Public Member Functions

 LinearizeSpeciesInterface (const GeomFields &geomFields, const T_Scalar RRConstant, const T_Scalar A_coeff, const T_Scalar B_coeff, const T_Scalar interfaceUnderRelax, const bool Anode, const bool Cathode, Field &varField, Field &mFElectricModel, Field &elecPotentialField)
 
void discretize (const Mesh &mesh, const Mesh &parentMesh, const Mesh &otherMesh, MultiFieldMatrix &mfmatrix, MultiField &xField, MultiField &rField)
 

Private Attributes

const GeomFields & _geomFields
 
Field & _varField
 
Field & _mFElectricModel
 
Field & _elecPotentialField
 
const T_Scalar _RRConstant
 
const T_Scalar _A_coeff
 
const T_Scalar _B_coeff
 
const T_Scalar _interfaceUnderRelax
 
const bool _Anode
 
const bool _Cathode
 

Detailed Description

template<class X, class Diag, class OffDiag>
class LinearizeSpeciesInterface< X, Diag, OffDiag >

Definition at line 24 of file LinearizeSpeciesInterface.h.

Member Typedef Documentation

template<class X, class Diag, class OffDiag>
typedef CRMatrix<Diag,OffDiag,X> LinearizeSpeciesInterface< X, Diag, OffDiag >::CCMatrix

Definition at line 28 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
typedef CCMatrix::DiagArray LinearizeSpeciesInterface< X, Diag, OffDiag >::DiagArray

Definition at line 29 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
typedef NumTypeTraits<X>::T_Scalar LinearizeSpeciesInterface< X, Diag, OffDiag >::T_Scalar

Definition at line 27 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
typedef Array<T_Scalar> LinearizeSpeciesInterface< X, Diag, OffDiag >::TArray

Definition at line 30 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
typedef Vector<T_Scalar,3> LinearizeSpeciesInterface< X, Diag, OffDiag >::VectorT3

Definition at line 32 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
typedef Array<VectorT3> LinearizeSpeciesInterface< X, Diag, OffDiag >::VectorT3Array

Definition at line 33 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
typedef Array<X> LinearizeSpeciesInterface< X, Diag, OffDiag >::XArray

Definition at line 31 of file LinearizeSpeciesInterface.h.

Constructor & Destructor Documentation

template<class X, class Diag, class OffDiag>
LinearizeSpeciesInterface< X, Diag, OffDiag >::LinearizeSpeciesInterface ( const GeomFields &  geomFields,
const T_Scalar  RRConstant,
const T_Scalar  A_coeff,
const T_Scalar  B_coeff,
const T_Scalar  interfaceUnderRelax,
const bool  Anode,
const bool  Cathode,
Field &  varField,
Field &  mFElectricModel,
Field &  elecPotentialField 
)
inline

Definition at line 36 of file LinearizeSpeciesInterface.h.

45  :
46  _geomFields(geomFields),
47  _varField(varField),
48  _mFElectricModel(mFElectricModel),
49  _elecPotentialField(elecPotentialField),
50  _RRConstant(RRConstant),
51  _A_coeff(A_coeff),
52  _B_coeff(B_coeff),
53  _interfaceUnderRelax(interfaceUnderRelax),
54  _Anode(Anode),
55  _Cathode(Cathode)
56  {}

Member Function Documentation

template<class X, class Diag, class OffDiag>
void LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize ( const Mesh &  mesh,
const Mesh &  parentMesh,
const Mesh &  otherMesh,
MultiFieldMatrix &  mfmatrix,
MultiField &  xField,
MultiField &  rField 
)
inline

Definition at line 59 of file LinearizeSpeciesInterface.h.

References LinearizeSpeciesInterface< X, Diag, OffDiag >::_Anode, LinearizeSpeciesInterface< X, Diag, OffDiag >::_Cathode, LinearizeSpeciesInterface< X, Diag, OffDiag >::_elecPotentialField, LinearizeSpeciesInterface< X, Diag, OffDiag >::_geomFields, LinearizeSpeciesInterface< X, Diag, OffDiag >::_interfaceUnderRelax, LinearizeSpeciesInterface< X, Diag, OffDiag >::_mFElectricModel, LinearizeSpeciesInterface< X, Diag, OffDiag >::_RRConstant, LinearizeSpeciesInterface< X, Diag, OffDiag >::_varField, GeomFields::areaMag, Mesh::getCellCells(), Mesh::getCells(), CRMatrix< T_Diag, T_OffDiag, X >::getCoeff(), StorageSite::getCount(), CRMatrix< T_Diag, T_OffDiag, X >::getDiag(), Mesh::getFaceCells(), MultiFieldMatrix::getMatrix(), Mesh::getOtherFaceGroupSite(), Mesh::getParentFaceGroupSite(), and StorageSite::getSelfCount().

Referenced by SpeciesModel< T >::Impl::linearize().

62  {
63  const StorageSite& cells = mesh.getCells();
64  const StorageSite& faces = mesh.getParentFaceGroupSite();
65 
66  // previous timestep shell mesh info
67  const XArray& mFElecModel =
68  dynamic_cast<const XArray&>(_mFElectricModel[cells]);
69 
70  // shell mesh info
71  const MultiField::ArrayIndex cVarIndex(&_varField,&cells);
72  CCMatrix& matrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndex,cVarIndex));
73  const XArray& xCell = dynamic_cast<const XArray&>(xField[cVarIndex]);
74  const CRConnectivity& cellCells = mesh.getCellCells();
75  XArray& rCell = dynamic_cast<XArray&>(rField[cVarIndex]);
76  DiagArray& diag = matrix.getDiag();
77 
78  // shell mesh electrical potential values
79  const XArray& ePotCell =
80  dynamic_cast<const XArray&>(_elecPotentialField[cells]);
81 
82 
83  // In the following, parent is assumed to be the electrolyte, and
84  // the other mesh is assumed to be electrode when implimenting
85  // the Butler-Volmer equations
86 
87  // parent mesh info
88  const CRConnectivity& parentFaceCells = parentMesh.getFaceCells(faces);
89  const StorageSite& parentCells = parentMesh.getCells();
90  const MultiField::ArrayIndex cVarIndexParent(&_varField,&parentCells);
91  XArray& rParentCell = dynamic_cast<XArray&>(rField[cVarIndexParent]);
92  CCMatrix& parentmatrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndexParent,cVarIndexParent));
93  DiagArray& parentdiag = parentmatrix.getDiag();
94  const TArray& faceAreaMag =
95  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
96 
97  // other mesh info
98  const StorageSite& otherFaces = mesh.getOtherFaceGroupSite();
99  const CRConnectivity& otherFaceCells = otherMesh.getFaceCells(otherFaces);
100  const StorageSite& otherCells = otherMesh.getCells();
101  const MultiField::ArrayIndex cVarIndexOther(&_varField,&otherCells);
102  XArray& rOtherCell = dynamic_cast<XArray&>(rField[cVarIndexOther]);
103  CCMatrix& othermatrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndexOther,cVarIndexOther));
104  DiagArray& otherdiag = othermatrix.getDiag();
105 
106  // set constants for entire shell
107  const T_Scalar F = 96485.0; // C/mol
108  const T_Scalar k = _RRConstant;
109  T_Scalar csMax = 26000.0;
110  if (_Anode){
111  csMax = 26390.0;}
112  if (_Cathode){
113  csMax = 22860.0;}
114  const T_Scalar alpha_a = 0.5;
115  const T_Scalar alpha_c = 0.5;
116  const T_Scalar R = 8.314; // J/mol/K
117  const T_Scalar Temp = 300.0; // K
118  const T_Scalar C_a = alpha_a*F/R/Temp;
119  const T_Scalar C_c = alpha_c*F/R/Temp;
120 
121  for (int f=0; f<faces.getCount(); f++)
122  {
123  //get parent mesh fluxes and coeffs
124  int c0p = parentFaceCells(f,0);
125  int c1p = parentFaceCells(f,1);
126  if (c1p < parentCells.getSelfCount())
127  {
128  // c0 is ghost cell and c1 is boundry cell, so swap cell numbers
129  // so that c1p refers to the ghost cell in the following
130  int temp = c0p;
131  c0p = c1p;
132  c1p = temp;
133  }
134 
135  const X parentFlux = rParentCell[c1p]; // inward shell flux on the left
136  const OffDiag dRC0dXC3 = parentmatrix.getCoeff(c1p, c0p);
137  const Diag dRC0dXC0 = parentdiag[c1p];
138 
139  //get other mesh fluxes and coeffs
140  int c0o = otherFaceCells(f,0);
141  int c1o = otherFaceCells(f,1);
142  if (c1o < otherCells.getSelfCount())
143  {
144  // c0 is ghost cell and c1 is boundry cell, so swap cell numbers
145  // so that c1o refers to the ghost cell in the following
146  int temp = c0o;
147  c0o = c1o;
148  c1o = temp;
149  }
150 
151  const X otherFlux = rOtherCell[c1o]; // inward shell flux on the right
152  const OffDiag dRC0dXC2 = othermatrix.getCoeff(c1o, c0o);
153  const OffDiag dRC0dXC1 = otherdiag[c1o];
154 
155 
156  //now put flux information from meshes into shell cells
157  const int c0 = f;
158  const int c1 = cellCells(f,0);
159  const int c2 = cellCells(f,1);
160  const int c3 = cellCells(f,2);
161 
162  const T_Scalar Area = faceAreaMag[c0];
163 
164  T_Scalar cs_star = xCell[c1];
165  T_Scalar ce_star = xCell[c0];
166 
167  // avoid nans
168  if (cs_star < 0.0){ cout << "ERROR: Cs < 0, Cs=" << cs_star << endl; cs_star = 0.0;}
169  if (ce_star < 0.0){ cout << "ERROR: Ce < 0, Ce=" << ce_star << endl; ce_star = 0.0;}
170  if (cs_star > csMax){ cout << "ERROR: Cs > CsMax, Cs=" << cs_star << endl; cs_star = csMax;}
171 
172  const T_Scalar SOC = mFElecModel[c1]/csMax;
173 
174  T_Scalar U_ref = 0.1; // V
175  if (_Anode){
176  U_ref = -0.16 + 1.32*exp(-3.0*SOC)+10.0*exp(-2000.0*SOC);
177  if (U_ref < 0.0)
178  {U_ref = 0.0; cout << "U_ref < 0" << endl;}
179  if (U_ref > 1.2)
180  {U_ref = 1.2; cout << "U_ref > 1.2" << endl;}
181  }
182  if (_Cathode){
183  U_ref = 4.19829 + 0.0565661*tanh(-14.5546*SOC + 8.60942) - 0.0275479*(1.0/pow((0.998432-SOC),0.492465) - 1.90111) - 0.157123*exp(-0.04738*pow(SOC,8.0)) + 0.810239*exp(-40.0*(SOC-0.133875));
184  if (U_ref < 0.0)
185  {U_ref = 0.0; cout << "U_ref < 0" << endl;}
186  if (U_ref > 5.0)
187  {U_ref = 5.0; cout << "U_ref > 5.0" << endl;}
188  }
189 
190  const T_Scalar eta_star = ePotCell[c1] - ePotCell[c0] - U_ref;
191 
192  //const T_Scalar eta_star = _A_coeff - _B_coeff - U_ref;
193  T_Scalar C_0 = exp(C_a*eta_star)-exp(-1.0*C_c*eta_star);
194 
195  const T_Scalar i_star = C_0*k*F*Area*pow(ce_star,alpha_c)*pow((csMax-cs_star),alpha_a)*pow(cs_star,alpha_c);
196 
197  // CURRENT SHOULD NOT BE ZERO
198  if (i_star == 0.0){cout << "WARNING: current = 0" << endl;}
199 
200  // calculate dC_0/dCS
201  T_Scalar dC_0dCS = 0.0;
202 
203  if (_Anode)
204  {
205  dC_0dCS = (C_a*exp(C_a*eta_star) + C_c*exp(-1*C_c*eta_star))*(-1.0)*(-20000.0*exp(-2000.0*SOC) - 3.96*exp(-3.0*SOC))*(1.0/csMax);
206  }
207  if (_Cathode)
208  {
209  dC_0dCS = (C_a*exp(C_a*eta_star) + C_c*exp(-1*C_c*eta_star))*(-1.0)*(-0.0135664/pow((0.998432-SOC),1.49247) - 0.823297/pow(cosh(8.60942-14.5546*SOC),2.0) + 0.0595559*exp(-0.04738*pow(SOC,8.0))*pow(SOC,7.0) - 6859.94*exp(-40.0*SOC))*(1.0/csMax);
210  }
211 
212  const T_Scalar dIdCS_star = i_star*(alpha_c/cs_star - alpha_a/(csMax-cs_star)+ dC_0dCS/C_0);
213 
214  const T_Scalar dIdCE_star = i_star*alpha_c/ce_star;
215 
216  // left(parent) shell cell - 3 neighbors
217  // flux balance
218  OffDiag& offdiagC0_C1 = matrix.getCoeff(c0, c1);
219  OffDiag& offdiagC0_C2 = matrix.getCoeff(c0, c2);
220  OffDiag& offdiagC0_C3 = matrix.getCoeff(c0, c3);
221 
222  rCell[c0] = F*otherFlux + F*parentFlux;
223  offdiagC0_C1 = F*dRC0dXC1;
224  offdiagC0_C3 = F*dRC0dXC3;
225  offdiagC0_C2 = F*dRC0dXC2;
226  diag[c0] = F*dRC0dXC0;
227 
228  // right(other) shell cell - 2 neighbors
229  // jump condition
230  OffDiag& offdiagC1_C0 = matrix.getCoeff(c1, c0);
231  OffDiag& offdiagC1_C2 = matrix.getCoeff(c1, c2);
232 
233  rCell[c1] = otherFlux*F - (i_star + dIdCS_star*(xCell[c1]-cs_star) + dIdCE_star*(xCell[c0]-ce_star));
234  offdiagC1_C0 = -1*dIdCE_star;
235  offdiagC1_C2 = F*dRC0dXC2;
236  //make sure diag is < 0
237  if (dIdCS_star > 0.0)
238  {
239  diag[c1] = F*dRC0dXC1 - dIdCS_star;
240  }
241  else
242  {
243  diag[c1] = F*dRC0dXC1;
244  }
245 
246  // set other coeffs to zero for right shell cell
247  OffDiag& offdiagC1_C3 = matrix.getCoeff(c1, c3);
248  offdiagC1_C3 = 0.0;
249 
250  // underrelax diagonal to help convergence
251  diag[c1] = 1/_interfaceUnderRelax*diag[c1];
252 
253  // some output of prevailing or starred values - useful for testing
254 
255  /*if (c0 == 0){
256  cout << "C_0: " << C_0 << endl;
257  cout << "C_s: " << cs_star << " C_e: " << ce_star << " i: " << i_star << " dIdCs: " << dIdCS_star << " dIdCe: " << dIdCE_star << endl;
258  cout <<"Diag: " << diag[c1] << " dRdXc2: " << offdiagC1_C2 << " dRdXc0: " << offdiagC1_C0 << endl;
259  cout << " " << endl;
260  }*/
261 
262  if (c0 == 0){
263  if (_Cathode)
264  {
265  //cout << "UrefAnode: " << U_ref << endl;
266  cout << "Cs0: " << cs_star << endl;
267  }
268  }
269 
270  }
271 
272  }
MultiFieldMatrix::getMatrix
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
Definition: MultiFieldMatrix.cpp:42
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
LinearizeSpeciesInterface::CCMatrix
CRMatrix< Diag, OffDiag, X > CCMatrix
Definition: LinearizeSpeciesInterface.h:28
Mesh::getParentFaceGroupSite
const StorageSite & getParentFaceGroupSite() const
Definition: Mesh.h:329
Mesh::getOtherFaceGroupSite
const StorageSite & getOtherFaceGroupSite() const
Definition: Mesh.h:332
MultiField::ArrayIndex
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
Mesh::getCellCells
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
LinearizeSpeciesInterface::T_Scalar
NumTypeTraits< X >::T_Scalar T_Scalar
Definition: LinearizeSpeciesInterface.h:27
Mesh::getCells
const StorageSite & getCells() const
Definition: Mesh.h:109
Mesh::getFaceCells
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
GeomFields::areaMag
Field areaMag
Definition: GeomFields.h:25

Member Data Documentation

template<class X, class Diag, class OffDiag>
const T_Scalar LinearizeSpeciesInterface< X, Diag, OffDiag >::_A_coeff
private

Definition at line 280 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
const bool LinearizeSpeciesInterface< X, Diag, OffDiag >::_Anode
private

Definition at line 283 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
const T_Scalar LinearizeSpeciesInterface< X, Diag, OffDiag >::_B_coeff
private

Definition at line 281 of file LinearizeSpeciesInterface.h.

template<class X, class Diag, class OffDiag>
const bool LinearizeSpeciesInterface< X, Diag, OffDiag >::_Cathode
private

Definition at line 284 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
Field& LinearizeSpeciesInterface< X, Diag, OffDiag >::_elecPotentialField
private

Definition at line 278 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
const GeomFields& LinearizeSpeciesInterface< X, Diag, OffDiag >::_geomFields
private

Definition at line 275 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
const T_Scalar LinearizeSpeciesInterface< X, Diag, OffDiag >::_interfaceUnderRelax
private

Definition at line 282 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
Field& LinearizeSpeciesInterface< X, Diag, OffDiag >::_mFElectricModel
private

Definition at line 277 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
const T_Scalar LinearizeSpeciesInterface< X, Diag, OffDiag >::_RRConstant
private

Definition at line 279 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

template<class X, class Diag, class OffDiag>
Field& LinearizeSpeciesInterface< X, Diag, OffDiag >::_varField
private

Definition at line 276 of file LinearizeSpeciesInterface.h.

Referenced by LinearizeSpeciesInterface< X, Diag, OffDiag >::discretize().

The documentation for this class was generated from the following file: