65 const XArray& xCell =
dynamic_cast<const XArray&
>(xField[cVarIndex]);
67 XArray& rCell =
dynamic_cast<XArray&
>(rField[cVarIndex]);
71 const XArray& eSpecConcCell =
83 XArray& rParentCell =
dynamic_cast<XArray&
>(rField[cVarIndexParent]);
85 DiagArray& parentdiag = parentmatrix.getDiag();
86 const TArray& faceAreaMag =
94 XArray& rOtherCell =
dynamic_cast<XArray&
>(rField[cVarIndexOther]);
96 DiagArray& otherdiag = othermatrix.getDiag();
111 const T_Scalar C_a = alpha_a*F/R/Temp;
112 const T_Scalar C_c = alpha_c*F/R/Temp;
115 for (
int f=0; f<faces.
getCount(); f++)
118 int c0p = parentFaceCells(f,0);
119 int c1p = parentFaceCells(f,1);
129 const X parentFlux = rParentCell[c1p];
130 const OffDiag dRC0dXC3 = parentmatrix.getCoeff(c1p, c0p);
131 const Diag dRC0dXC0 = parentdiag[c1p];
134 int c0o = otherFaceCells(f,0);
135 int c1o = otherFaceCells(f,1);
145 const X otherFlux = rOtherCell[c1o];
146 const OffDiag dRC0dXC2 = othermatrix.getCoeff(c1o, c0o);
147 const OffDiag dRC0dXC1 = otherdiag[c1o];
152 const int c1 = cellCells(f,0);
153 const int c2 = cellCells(f,1);
154 const int c3 = cellCells(f,2);
156 const T_Scalar Area = faceAreaMag[c0];
157 T_Scalar Ce_star = eSpecConcCell[c0];
158 T_Scalar Cs_star = eSpecConcCell[c1];
165 cout <<
"ERROR: Ce_star < 0, Ce_star=" << Ce_star << endl; Ce_star = 0.0;
169 cout <<
"ERROR: Cs_star < 0, Cs_star=" << Cs_star << endl; Cs_star = 0.0;
171 if (Cs_star > csMax){ Cs_star = 0.9*csMax; cout <<
"ERROR: Cs > CsMax" << endl;}
177 U_ref = -0.16 + 1.32*exp(-3.0*SOC)+10.0*exp(-2000.0*SOC);}
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));}
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);
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;
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));
195 OffDiag& offdiagC0_C1 = matrix.getCoeff(c0, c1);
196 OffDiag& offdiagC0_C2 = matrix.getCoeff(c0, c2);
197 OffDiag& offdiagC0_C3 = matrix.getCoeff(c0, c3);
199 rCell[c0] = otherFlux + parentFlux;
200 offdiagC0_C1 = dRC0dXC1;
201 offdiagC0_C3 = dRC0dXC3;
202 offdiagC0_C2 = dRC0dXC2;
207 OffDiag& offdiagC1_C0 = matrix.getCoeff(c1, c0);
208 OffDiag& offdiagC1_C2 = matrix.getCoeff(c1, c2);
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;
218 cout <<
"Warning: Diag > 0" << endl;
222 OffDiag& offdiagC1_C3 = matrix.getCoeff(c1, c3);
Matrix & getMatrix(const Index &rowIndex, const Index &colIndex)
CCMatrix::DiagArray DiagArray
const StorageSite & getParentFaceGroupSite() const
const GeomFields & _geomFields
CRMatrix< Diag, OffDiag, X > CCMatrix
Field & _speciesConcentrationField
const StorageSite & getOtherFaceGroupSite() const
const T_Scalar _RRConstant
pair< const Field *, const StorageSite * > ArrayIndex
const CRConnectivity & getCellCells() const
const StorageSite & getCells() const
const CRConnectivity & getFaceCells(const StorageSite &site) const
NumTypeTraits< X >::T_Scalar T_Scalar