Memosa-FVM  0.2
LinearizePotentialInterface< X, Diag, OffDiag > Class Template Reference

#include <LinearizePotentialInterface.h>

Collaboration diagram for LinearizePotentialInterface< X, Diag, OffDiag >:

Public Types

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

Public Member Functions

 LinearizePotentialInterface (const GeomFields &geomFields, Field &varField, Field &speciesConcentrationField, const T_Scalar RRConstant, const T_Scalar A_coeff, const T_Scalar B_coeff, const bool Anode, const bool Cathode)
 
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_speciesConcentrationField
 
const T_Scalar _RRConstant
 
const T_Scalar _A_coeff
 
const T_Scalar _B_coeff
 
const bool _Anode
 
const bool _Cathode
 

Detailed Description

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

Definition at line 24 of file LinearizePotentialInterface.h.

Member Typedef Documentation

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

Definition at line 28 of file LinearizePotentialInterface.h.

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

Definition at line 29 of file LinearizePotentialInterface.h.

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

Definition at line 27 of file LinearizePotentialInterface.h.

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

Definition at line 30 of file LinearizePotentialInterface.h.

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

Definition at line 32 of file LinearizePotentialInterface.h.

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

Definition at line 33 of file LinearizePotentialInterface.h.

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

Definition at line 31 of file LinearizePotentialInterface.h.

Constructor & Destructor Documentation

template<class X, class Diag, class OffDiag>
LinearizePotentialInterface< X, Diag, OffDiag >::LinearizePotentialInterface ( const GeomFields geomFields,
Field varField,
Field speciesConcentrationField,
const T_Scalar  RRConstant,
const T_Scalar  A_coeff,
const T_Scalar  B_coeff,
const bool  Anode,
const bool  Cathode 
)
inline

Definition at line 36 of file LinearizePotentialInterface.h.

43  :
44  _geomFields(geomFields),
45  _varField(varField),
46  _speciesConcentrationField(speciesConcentrationField),
47  _RRConstant(RRConstant),
48  _A_coeff(A_coeff),
49  _B_coeff(B_coeff),
50  _Anode(Anode),
51  _Cathode(Cathode)
52  {}

Member Function Documentation

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

Definition at line 55 of file LinearizePotentialInterface.h.

References LinearizePotentialInterface< X, Diag, OffDiag >::_Anode, LinearizePotentialInterface< X, Diag, OffDiag >::_Cathode, LinearizePotentialInterface< X, Diag, OffDiag >::_geomFields, LinearizePotentialInterface< X, Diag, OffDiag >::_RRConstant, LinearizePotentialInterface< X, Diag, OffDiag >::_speciesConcentrationField, LinearizePotentialInterface< 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 ElectricModel< T >::Impl::linearizeElectroStatics().

58  {
59  const StorageSite& cells = mesh.getCells();
60  const StorageSite& faces = mesh.getParentFaceGroupSite();
61 
62  // shell mesh info
63  const MultiField::ArrayIndex cVarIndex(&_varField,&cells);
64  CCMatrix& matrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndex,cVarIndex));
65  const XArray& xCell = dynamic_cast<const XArray&>(xField[cVarIndex]);
66  const CRConnectivity& cellCells = mesh.getCellCells();
67  XArray& rCell = dynamic_cast<XArray&>(rField[cVarIndex]);
68  DiagArray& diag = matrix.getDiag();
69 
70  //shell mesh species concentration values
71  const XArray& eSpecConcCell =
72  dynamic_cast<const XArray&>(_speciesConcentrationField[cells]);
73 
74 
75  // In the following, parent is assumed to be the electrolyte, and
76  // the other mesh is assumed to be electrode when implimenting
77  // the Butler-Volmer equations
78 
79  // parent mesh info
80  const CRConnectivity& parentFaceCells = parentMesh.getFaceCells(faces);
81  const StorageSite& parentCells = parentMesh.getCells();
82  const MultiField::ArrayIndex cVarIndexParent(&_varField,&parentCells);
83  XArray& rParentCell = dynamic_cast<XArray&>(rField[cVarIndexParent]);
84  CCMatrix& parentmatrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndexParent,cVarIndexParent));
85  DiagArray& parentdiag = parentmatrix.getDiag();
86  const TArray& faceAreaMag =
87  dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);
88 
89  // other mesh info
90  const StorageSite& otherFaces = mesh.getOtherFaceGroupSite();
91  const CRConnectivity& otherFaceCells = otherMesh.getFaceCells(otherFaces);
92  const StorageSite& otherCells = otherMesh.getCells();
93  const MultiField::ArrayIndex cVarIndexOther(&_varField,&otherCells);
94  XArray& rOtherCell = dynamic_cast<XArray&>(rField[cVarIndexOther]);
95  CCMatrix& othermatrix = dynamic_cast<CCMatrix&>(mfmatrix.getMatrix(cVarIndexOther,cVarIndexOther));
96  DiagArray& otherdiag = othermatrix.getDiag();
97 
98  // set constants for entire shell
99  const T_Scalar F = 96485.0; // C/mol
100  const T_Scalar k = _RRConstant;
101  T_Scalar csMax = 26000.0;// mol/m^3
102  if (_Anode){
103  csMax = 26390.0;}
104  if (_Cathode){
105  csMax = 22860.0;}
106 
107  const T_Scalar alpha_a = 0.5;
108  const T_Scalar alpha_c = 0.5;
109  const T_Scalar R = 8.314; // J/mol/K
110  const T_Scalar Temp = 300.0; // K
111  const T_Scalar C_a = alpha_a*F/R/Temp;
112  const T_Scalar C_c = alpha_c*F/R/Temp;
113 
114 
115  for (int f=0; f<faces.getCount(); f++)
116  {
117  //get parent mesh fluxes and coeffs
118  int c0p = parentFaceCells(f,0);
119  int c1p = parentFaceCells(f,1);
120  if (c1p < parentCells.getSelfCount())
121  {
122  // c0 is ghost cell and c1 is boundry cell, so swap cell numbers
123  // so that c1p refers to the ghost cell in the following
124  int temp = c0p;
125  c0p = c1p;
126  c1p = temp;
127  }
128 
129  const X parentFlux = rParentCell[c1p]; // inward shell flux on the left
130  const OffDiag dRC0dXC3 = parentmatrix.getCoeff(c1p, c0p);
131  const Diag dRC0dXC0 = parentdiag[c1p];
132 
133  //get other mesh fluxes and coeffs
134  int c0o = otherFaceCells(f,0);
135  int c1o = otherFaceCells(f,1);
136  if (c1o < otherCells.getSelfCount())
137  {
138  // c0 is ghost cell and c1 is boundry cell, so swap cell numbers
139  // so that c1o refers to the ghost cell in the following
140  int temp = c0o;
141  c0o = c1o;
142  c1o = temp;
143  }
144 
145  const X otherFlux = rOtherCell[c1o]; // inward shell flux on the right
146  const OffDiag dRC0dXC2 = othermatrix.getCoeff(c1o, c0o);
147  const OffDiag dRC0dXC1 = otherdiag[c1o];
148 
149 
150  //now put flux information from meshes into shell cells
151  const int c0 = f;
152  const int c1 = cellCells(f,0);
153  const int c2 = cellCells(f,1);
154  const int c3 = cellCells(f,2);
155 
156  const T_Scalar Area = faceAreaMag[c0];
157  T_Scalar Ce_star = eSpecConcCell[c0];
158  T_Scalar Cs_star = eSpecConcCell[c1];
159 
160 
161 
162  //avoid nans during iterations
163  if (Ce_star < 0)
164  {
165  cout << "ERROR: Ce_star < 0, Ce_star=" << Ce_star << endl; Ce_star = 0.0;
166  }
167  if (Cs_star < 0)
168  {
169  cout << "ERROR: Cs_star < 0, Cs_star=" << Cs_star << endl; Cs_star = 0.0;
170  }
171  if (Cs_star > csMax){ Cs_star = 0.9*csMax; cout << "ERROR: Cs > CsMax" << endl;}
172 
173  const T_Scalar SOC = Cs_star/csMax;
174 
175  T_Scalar U_ref = 0.1; // V
176  if (_Anode){
177  U_ref = -0.16 + 1.32*exp(-3.0*SOC)+10.0*exp(-2000.0*SOC);}
178  if (_Cathode){
179  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));}
180 
181  const T_Scalar i0_star = k*F*Area*pow(Ce_star,alpha_c)*pow((csMax-Cs_star),alpha_a)*pow(Cs_star,alpha_c);
182 
183  const T_Scalar Phis_star = xCell[c1];
184  const T_Scalar Phie_star = xCell[c0];
185  T_Scalar eta_star = Phis_star-Phie_star-U_ref;
186 
187  //cout << eta_star << endl;
188 
189  const T_Scalar i_star = i0_star*(exp(C_a*eta_star)-exp(-1*C_c*eta_star));
190  const T_Scalar dIdPhiS_star = i0_star*(C_a*exp(C_a*eta_star)+C_c*exp(-1*C_c*eta_star));
191  const T_Scalar dIdPhiE_star = -1*i0_star*(C_a*exp(C_a*eta_star)+C_c*exp(-1*C_c*eta_star));
192 
193  // left(parent) shell cell - 3 neighbors
194  // flux balance
195  OffDiag& offdiagC0_C1 = matrix.getCoeff(c0, c1);
196  OffDiag& offdiagC0_C2 = matrix.getCoeff(c0, c2);
197  OffDiag& offdiagC0_C3 = matrix.getCoeff(c0, c3);
198 
199  rCell[c0] = otherFlux + parentFlux;
200  offdiagC0_C1 = dRC0dXC1;
201  offdiagC0_C3 = dRC0dXC3;
202  offdiagC0_C2 = dRC0dXC2;
203  diag[c0] = dRC0dXC0;
204 
205  // right(other) shell cell - 2 neighbors
206  // BV(jump) condition
207  OffDiag& offdiagC1_C0 = matrix.getCoeff(c1, c0);
208  OffDiag& offdiagC1_C2 = matrix.getCoeff(c1, c2);
209 
210  rCell[c1] = otherFlux - (i_star + dIdPhiS_star*(xCell[c1]-Phis_star) + dIdPhiE_star*(xCell[c0]-Phie_star));
211  offdiagC1_C0 = -1*dIdPhiE_star;
212  offdiagC1_C2 = dRC0dXC2;
213  diag[c1] = dRC0dXC1 - dIdPhiS_star;
214 
215  //make sure diag is < 0
216  if (diag[c1] > 0.0)
217  {
218  cout << "Warning: Diag > 0" << endl;
219  }
220 
221  // set other coeffs to zero for right shell cell
222  OffDiag& offdiagC1_C3 = matrix.getCoeff(c1, c3);
223  offdiagC1_C3 = 0.0;
224 
225  // some output of prevailing values once per shell - useful for testing
226 
227  /*if (c0 == 0){
228  cout << "i0_star: " << i0_star << endl;
229  cout << "Phi_s: " << Phis_star << " Phi_e: " << Phie_star << " i: " << i_star << " Flux: " << otherFlux << " dIdPhis: " << dIdPhiS_star << " dIdPhie: " << dIdPhiE_star << endl;
230  cout <<"Diag: " << diag[c1] << " dRdXc2: " << offdiagC1_C2 << " dRdXc0: " << offdiagC1_C0 << endl;
231  cout << " " << endl;
232  }*/
233 
234  //cout << "ParentFlux: " << parentFlux << " otherFlux: " << otherFlux << endl;
235  }
236 
237  }
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
int getSelfCount() const
Definition: StorageSite.h:40
const StorageSite & getParentFaceGroupSite() const
Definition: Mesh.h:329
CRMatrix< Diag, OffDiag, X > CCMatrix
const StorageSite & getOtherFaceGroupSite() const
Definition: Mesh.h:332
pair< const Field *, const StorageSite * > ArrayIndex
Definition: MultiField.h:21
const CRConnectivity & getCellCells() const
Definition: Mesh.cpp:480
const StorageSite & getCells() const
Definition: Mesh.h:109
const CRConnectivity & getFaceCells(const StorageSite &site) const
Definition: Mesh.cpp:388
int getCount() const
Definition: StorageSite.h:39
NumTypeTraits< X >::T_Scalar T_Scalar
Field areaMag
Definition: GeomFields.h:25

Member Data Documentation

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

Definition at line 244 of file LinearizePotentialInterface.h.

template<class X, class Diag, class OffDiag>
const bool LinearizePotentialInterface< X, Diag, OffDiag >::_Anode
private
template<class X, class Diag, class OffDiag>
const T_Scalar LinearizePotentialInterface< X, Diag, OffDiag >::_B_coeff
private

Definition at line 245 of file LinearizePotentialInterface.h.

template<class X, class Diag, class OffDiag>
const bool LinearizePotentialInterface< X, Diag, OffDiag >::_Cathode
private
template<class X, class Diag, class OffDiag>
const GeomFields& LinearizePotentialInterface< X, Diag, OffDiag >::_geomFields
private
template<class X, class Diag, class OffDiag>
const T_Scalar LinearizePotentialInterface< X, Diag, OffDiag >::_RRConstant
private
template<class X, class Diag, class OffDiag>
Field& LinearizePotentialInterface< X, Diag, OffDiag >::_speciesConcentrationField
private
template<class X, class Diag, class OffDiag>
Field& LinearizePotentialInterface< X, Diag, OffDiag >::_varField
private

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