00001 #include "facefluxcompositional.h"
00002 #include "numericmethods.h"
00003
00004
00005 FaceFluxCompositional::FaceFluxCompositional(OrthoMesh &mesh,FlashCompositional &flash,const VecDouble &K,double v1,double v2,double sr1,double sr2,double grav)
00006 :FaceFluxFunction(3),m_flash(flash),m_mob(v1,v2,sr1,sr2),m_gravMob(v1,v2,sr1,sr2),m_Dmob(v1,v2,sr1,sr2),m_DgravMob(v1,v2,sr1,sr2),m_pP(NULL),m_K(K),m_grav(grav)
00007 {
00008 assert(flash.numPhases()==2);
00009 m_vNormalVel.reinit(mesh.numFaces());
00010 assert(m_grav>=0.);
00011
00012 m_densities1.reinit(2);
00013 m_densities2.reinit(2);
00014
00015 }
00016
00017
00018
00019 void FaceFluxCompositional::fluxAtFace(VecDouble &vFlux,OrthoMesh::Face_It &faceIt,int face,int cell1,int cell2,const VecDouble &Q1,const VecDouble &Q2)
00020 {
00021 assert(vFlux.size());
00022 assert((unsigned) face < m_vNormalVel.size());
00023 static FlashData data1(2,3),data2(2,3);
00024 static VecDouble phasesVol(2);
00025 static VecDouble vFlux1(3);
00026 static VecDouble vFlux2(3);
00027 const VecDouble &P = getPressure();
00028 double K1 = getK()(cell1);
00029 double K2 = getK()(cell2);
00030 double vel = m_vNormalVel(face);
00031
00032 m_flash.flash(P(cell1),Q1,data1);
00033 m_flash.flash(P(cell2),Q2,data2);
00034
00035 m_flash.getPhasesVolume(P(cell1),data1,phasesVol);
00036 m_Sw1=phasesVol(0)/(phasesVol(0)+phasesVol(1));
00037 data1.getPhasesDensities(m_densities1,m_flash.getComponentsMolarMass(),phasesVol);
00038
00039 double mobW=m_mob(m_Sw1);
00040 double gravFlux=K1*m_gravMob(m_Sw1)*(m_densities1(1)-m_densities1(0))*m_grav;
00041 double vel0 = mobW*vel + gravFlux;
00042 double vel1 = (1.0-mobW)*vel - gravFlux;
00043 vFlux(0)=data1.getMoles(0,0)/phasesVol(0)*vel0 + data1.getMoles(1,0)/phasesVol(1)*vel1;
00044 vFlux(1)=data1.getMoles(0,1)/phasesVol(0)*vel0 + data1.getMoles(1,1)/phasesVol(1)*vel1;
00045 vFlux(2)=data1.getMoles(0,2)/phasesVol(0)*vel0 + 0;
00046
00047
00048
00049 m_flash.getPhasesVolume(P(cell2),data2,phasesVol);
00050 m_Sw2=phasesVol(0)/(phasesVol(0)+phasesVol(1));
00051 data2.getPhasesDensities(m_densities2,m_flash.getComponentsMolarMass(),phasesVol);
00052
00053
00054 mobW=m_mob(m_Sw2);
00055 gravFlux = K2*m_gravMob(m_Sw2)*(m_densities2(1)-m_densities2(0))*m_grav;
00056 vel0=mobW*vel + gravFlux;
00057 vel0=(1-mobW)*vel - gravFlux;
00058
00059 vFlux(0)+=data2.getMoles(0,0)/phasesVol(0)*vel0 + data2.getMoles(1,0)/phasesVol(1)*vel1;
00060 vFlux(1)+=data2.getMoles(0,1)/phasesVol(0)*vel0 + data2.getMoles(1,1)/phasesVol(1)*vel1;
00061 vFlux(2)+=data2.getMoles(0,2)/phasesVol(0)*vel0 + 0;
00062
00063 assert(data2.getMoles(1,2) == 0.0);
00064 assert(data1.getMoles(1,2) == 0.0);
00065
00066 vFlux*=0.5;
00067
00068
00069
00070 m_cell1=cell1;
00071 m_cell2=cell2;
00072 }
00073
00074
00075 void FaceFluxCompositional::updateDynamicData(DynamicBase &dynMod)
00076 {
00077 dynMod.getNormalVelocityAtFaces(m_vNormalVel);
00078 m_pP=&dynMod.getPressureAtCells();
00079 }
00080
00081
00082
00083 void FaceFluxCompositional::maxLocalCharVelocity(VecDouble &vout, OrthoMesh::Face_It &faceIt,int face,int cell1,int cell2,const VecDouble &Q1, const VecDouble &Q2)
00084 {
00085 assert(m_cell1==cell1);
00086 assert(m_cell2==cell2);
00087 if (m_Sw1 > m_Sw2)
00088 std::swap(m_Sw1,m_Sw2);
00089 double grav1 = getK()(cell1)*fabs(m_densities1(0)-m_densities1(1));
00090 double grav2 = getK()(cell2)*fabs(m_densities2(0)-m_densities2(1));
00091 vout=
00092 m_vNormalVel(face)*m_Dmob.getMaxNorm(m_Sw1,m_Sw2) +
00093 m_DgravMob.getMaxNorm(m_Sw1,m_Sw2)*std::max(grav1,grav2);
00094
00095
00096 }
00097
00098
00099 void FaceFluxCompositional::transformToComponentMoles(VecDouble &v1)
00100 {
00101 assert(v1.size()==3);
00102 return;
00103 }
00104
00105
00106 void FaceFluxCompositional::maxGlobalCharVelocity(double vel[3])
00107 {
00108 double maxVel=NumericMethods::vectorMaxAbsValue(m_vNormalVel);
00109 vel[X]= maxVel*m_Dmob.getMaxNorm(0,1);
00110 vel[Y]=0;
00111 vel[Z]=0;
00112 return;
00113 }
