00001 #include "facefluxco2brine.h"
00002 #include "numericmethods.h"
00003 #include <fstream>
00004
00005 FaceFluxCO2Brine::FaceFluxCO2Brine(OrthoMesh &mesh,FlashCompositional &flash,const VecDouble &K,GeneralFunctionInterface &fMob,GeneralFunctionInterface &DfMob,GeneralFunctionInterface &fGrav,GeneralFunctionInterface &DfGrav,double grav)
00006 :FaceFluxFunction(3),m_flash(flash),
00007 m_mob (dynamic_cast<Function1D&>(fMob)),
00008 m_gravMob(dynamic_cast<Function1D&>(fGrav)),
00009 m_Dmob (dynamic_cast<Function1D&>(DfMob)),
00010 m_DgravMob(dynamic_cast<Function1D&>(DfGrav)),
00011 m_pP(NULL),m_K(K),m_grav(grav),m_mesh(mesh)
00012 {
00013 assert(flash.numPhases()==2);
00014 m_vNormalVel.reinit(mesh.numFaces());
00015 assert(m_grav>=0.);
00016
00017 m_densities1.reinit(2);
00018 m_densities2.reinit(2);
00019 std::ofstream file("mobGrav.dat");
00020 for (double x=0;x<=1;x+=0.01)
00021 {
00022 file << x << " " << m_gravMob(x) << std::endl;
00023 }
00024
00025 }
00026
00027 using NumericMethods::_div;
00028
00029 void FaceFluxCO2Brine::fluxAtFace(VecDouble &vFlux, const FaceInfo &face,const VecDouble &Q1,const VecDouble &Q2)
00030 {
00031 assert(vFlux.size());
00032 assert((unsigned) face.index < m_vNormalVel.size());
00033 static FlashData data_1(2,3,NULL),data_2(2,3,NULL);
00034 static FlashData data_3(2,3),data_4(2,3);
00035 static FlashData data1(2,3,NULL),data2(2,3,NULL);
00036 static VecDouble phasesVol1(2);
00037 static VecDouble phasesVol2(2);
00038 static VecDouble vFlux1(3);
00039 static VecDouble vFlux2(3);
00040 const VecDouble &P = getPressure();
00041 double K1 = getK()(face.cell1);
00042 double K2 = getK()(face.cell2);
00043 double vel = m_vNormalVel(face.index);
00044 double gravFlux1;
00045 double gravFlux2;
00046 double P1=P(face.cell1);
00047 double P2=P(face.cell2);
00048 if (!face.at_boundary)
00049 {
00050 m_flash.flash(face.cell1,data_1);
00051 m_flash.flash(face.cell2,data_2);
00052 data1.setData(data_1);
00053 data2.setData(data_2);
00054 }
00055 else
00056 {
00057 if (Q1 == Q2)
00058 {
00059 m_flash.flash(face.cell1,data_1);
00060 m_flash.flash(face.cell2,data_2);
00061 data1.setData(data_1);
00062 data2.setData(data_2);
00063 }
00064 else
00065 {
00066 m_flash.flash(P1,Q1,data_3);
00067 m_flash.flash(P2,Q2,data_4);
00068 data1.setData(data_3);
00069 data2.setData(data_4);
00070 }
00071 }
00072 m_flash.getPhasesVolume(P1,data1,phasesVol1);
00073 m_Sw1=phasesVol1(0);
00074 data1.getPhasesDensities(m_densities1,m_flash.getComponentsMolarMass(),phasesVol1);
00075
00076
00077 m_flash.getPhasesVolume(P2,data2,phasesVol2);
00078 m_Sw2=phasesVol2(0);
00079 data2.getPhasesDensities(m_densities2,m_flash.getComponentsMolarMass(),phasesVol2);
00080
00081
00082 if (face.normal == OrthoMesh::NORMAL_Z && !(face.at_boundary))
00083 {
00084 gravFlux1=K1*m_gravMob(m_Sw1)*(m_densities1(1)-m_densities1(0))*m_grav;
00085 gravFlux2=K2*m_gravMob(m_Sw2)*(m_densities2(1)-m_densities2(0))*m_grav;
00086
00087 }
00088 else
00089 {
00090 gravFlux1=gravFlux2=0.0;
00091 }
00092 double vel0 = m_mob(m_Sw1)*vel + gravFlux1;
00093 double vel1 = vel-vel0;
00094 vFlux(0)=_div(data1.getMoles(0,0),phasesVol1(0))*vel0 + _div(data1.getMoles(1,0),phasesVol1(1))*vel1;
00095 vFlux(1)=_div(data1.getMoles(0,1),phasesVol1(0))*vel0 + _div(data1.getMoles(1,1),phasesVol1(1))*vel1;
00096 vFlux(2)=_div(data1.getMoles(0,2),phasesVol1(0))*vel0 + 0;
00097
00098
00099
00100 vel0=m_mob(m_Sw2)*vel + gravFlux2;
00101 vel1=vel-vel0;
00102
00103 vFlux(0)+=_div(data2.getMoles(0,0),phasesVol2(0))*vel0 + _div(data2.getMoles(1,0),phasesVol2(1))*vel1;
00104 vFlux(1)+=_div(data2.getMoles(0,1),phasesVol2(0))*vel0 + _div(data2.getMoles(1,1),phasesVol2(1))*vel1;
00105 vFlux(2)+=_div(data2.getMoles(0,2),phasesVol2(0))*vel0 + 0;
00106
00107 assert(data2.getMoles(1,2) == 0.0);
00108 assert(data1.getMoles(1,2) == 0.0);
00109
00110 vFlux*=0.5;
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 m_cell1=face.cell1;
00126 m_cell2=face.cell2;
00127
00128 }
00129
00130
00131
00132
00133
00134
00135
00136 void FaceFluxCO2Brine::updateDynamicData(DynamicBase &dynMod)
00137 {
00138 dynMod.getNormalVelocityAtFaces(m_vNormalVel);
00139 m_pP=&dynMod.getPressureAtCells();
00140 }
00141
00142
00143 double FaceFluxCO2Brine::maxLocalCharVelocity(const FaceInfo &face,const VecDouble &Q1, const VecDouble &Q2)
00144 {
00145 assert(m_cell1==face.cell1);
00146 assert(m_cell2==face.cell2);
00147 if (m_Sw1 > m_Sw2)
00148 std::swap(m_Sw1,m_Sw2);
00149 double grav1 = getK()(face.cell1)*fabs(m_densities1(0)-m_densities1(1))*m_grav;
00150 double grav2 = getK()(face.cell2)*fabs(m_densities2(0)-m_densities2(1))*m_grav;
00151
00152 double mobMaxNormGrad=m_Dmob.getMaxNorm(m_Sw1,m_Sw2);
00153 if (mobMaxNormGrad == 0)
00154 mobMaxNormGrad=1;
00155 if (face.normal == OrthoMesh::NORMAL_Z)
00156 return fabs(m_vNormalVel(face.index)*mobMaxNormGrad) + fabs(m_DgravMob.getMaxNorm(m_Sw1,m_Sw2)*std::max(grav1,grav2));
00157 else
00158 return fabs(m_vNormalVel(face.index)*mobMaxNormGrad);
00159
00160 }
00161
00162
00163
00164
00165
00166
00167
00168 void FaceFluxCO2Brine::transformToComponentMoles(VecDouble &v1)
00169 {
00170 assert(v1.size()==3);
00171 return;
00172 }
00173
00174
00175 void FaceFluxCO2Brine::maxGlobalCharVelocity(double vel[3])
00176 {
00177 static bool bFirstTime = true;
00178 static double gMaxMobGrad;
00179 static double gMaxGravTerm;
00180 static VecDouble Vdt(3);
00181 if (bFirstTime)
00182 {
00183 double minP=NumericMethods::vectorMinValue(getPressure());
00184 double maxK=NumericMethods::vectorMaxValue(getK());
00185 bFirstTime=false;
00186 VecDouble totalComponents(3),phasesVolume(2);
00187 VecDouble vDensities(2);
00188 FlashData data(m_flash.numPhases(),m_flash.numComponents());
00189 totalComponents(0)=1;
00190 totalComponents(1)=1;
00191 totalComponents(2)=0;
00192 m_flash.flash(minP,totalComponents,data);
00193 m_flash.getPhasesVolume(minP,data,phasesVolume);
00194 data.getPhasesDensities(vDensities,m_flash.getComponentsMolarMass(),phasesVolume);
00195 assert(vDensities(0) != 0.0 && vDensities(1) != 0.0);
00196
00197 gMaxMobGrad=m_Dmob.getMaxNorm(0,1);
00198 gMaxGravTerm=fabs(maxK*(vDensities(0)-vDensities(1))*m_DgravMob.getMaxNorm(0,1)*m_grav);
00199
00200 printf("Max char velocity calculation\n");
00201 printf("\tDensities %g %g\n",vDensities(0),vDensities(1));
00202 printf("\tMax K %g\n",maxK);
00203 printf("\tMax Mobility Grad %g\n",gMaxMobGrad);
00204 printf("\tMax Grad Grav Term %g\n",gMaxGravTerm);
00205 printf("\tMax Grad Mobm %g\n",m_DgravMob.getMaxNorm(0,1));
00206
00207
00208 }
00209
00210
00211 m_mesh.getMaxAbsVelocitiesInFacesByComponents(Vdt,m_vNormalVel);
00212
00213 vel[X]= Vdt(X)*gMaxMobGrad;
00214 vel[Y]= Vdt(Y)*gMaxMobGrad;
00215 vel[Z]= Vdt(Z)*gMaxMobGrad + gMaxGravTerm;
00216
00217 printf("Char Vel %g %g %g\n",vel[X],vel[Y],vel[Z]);
00218 return;
00219 }
00220
00221
00222
00223
00224
00225
