00001 #include "facefluxsimpleblackoil.h"
00002 #include "numericmethods.h"
00003
00004 FaceFluxSimpleBlackOil::FaceFluxSimpleBlackOil(OrthoMesh &mesh,FlashCompositional &flash) :FaceFluxFunction(2),m_mesh(mesh),m_flash(flash),m_mob(1,1,0,0),m_Dmob(1,1,0,0),m_pP(NULL),m_vViscosities1(2),m_vViscosities2(2)
00005 {
00006 assert(flash.numPhases()==2);
00007 assert(flash.numComponents() == 2);
00008 m_vNormalVel.reinit(mesh.numFaces());
00009 }
00010
00011
00012
00013
00014 void FaceFluxSimpleBlackOil::updateDynamicData(DynamicBase &dynMod)
00015 {
00016 dynMod.getNormalVelocityAtFaces(m_vNormalVel);
00017 m_pP=&(dynMod.getPressureAtCells());
00018 }
00019
00020
00021
00022
00023
00024 void FaceFluxSimpleBlackOil::maxGlobalCharVelocity(double vel[3])
00025 {
00026 double maxVel=NumericMethods::vectorMaxAbsValue(m_vNormalVel);
00027 double minP=NumericMethods::vectorMinAbsValue(getPressure());
00028 static FlashData data(2,2);
00029 data.allocateOwnMemory();
00030 m_flash.getPhasesViscosities(minP,data,m_vViscosities1);
00031 m_Dmob.setParameters(m_vViscosities1);
00032 m_Dmob.updateInflectionPoint();
00033 vel[X]= maxVel*m_Dmob.getMaxNorm(0,1);
00034 vel[Y]=0;
00035 vel[Z]=0;
00036 return;
00037
00038 }
00039
00040
00041 void FaceFluxSimpleBlackOil::transformToFlash(VecDouble &vec)
00042 {
00043 return;
00044 }
00045
00046
00047 void FaceFluxSimpleBlackOil::fluxAtFace(VecDouble &vFlux, const FaceInfo &face,const VecDouble &Q1,const VecDouble &Q2)
00048 {
00049 assert(vFlux.size()==2);
00050 assert((unsigned) face.index < m_vNormalVel.size());
00051 static FlashData data1(2,2),data2(2,2);
00052 static VecDouble phasesVolume(2);
00053 static VecDouble vFlux1(2);
00054 static VecDouble vFlux2(2);
00055 const VecDouble &P = getPressure();
00056 m_bUnsaturated=false;
00057
00058 m_flash.flash(P(face.cell1),Q1,data1);
00059 m_flash.flash(P(face.cell2),Q2,data2);
00060
00061 m_flash.getPhasesViscosities(P(face.cell1),data1,m_vViscosities1);
00062 m_flash.getPhasesViscosities(P(face.cell2),data2,m_vViscosities2);
00063
00064 m_flash.getPhasesVolume(P(face.cell1),data1,phasesVolume);
00065 m_Sw1=phasesVolume(0);
00066
00067
00068 m_mob.setParameters(m_vViscosities1);
00069 double mobW=m_mob(m_Sw1);
00070 vFlux(0)=0.0;
00071 vFlux(1)=0.0;
00072
00073 if (phasesVolume(0) != 0.0)
00074 {
00075 m_cl11=data1.getMoles(0,0)/phasesVolume(0);
00076 m_cl21=data1.getMoles(0,1)/phasesVolume(0);
00077 vFlux(0)=data1.getMoles(0,0)/phasesVolume(0)*mobW;
00078 vFlux(1)=data1.getMoles(0,1)/phasesVolume(0)*mobW;
00079 }
00080 else
00081 m_cl11=m_cl21=0;
00082
00083 if (phasesVolume(1) != 0.0)
00084 {
00085 m_cl12=data1.getMoles(1,0)/phasesVolume(1);
00086 m_cl22=data1.getMoles(1,1)/phasesVolume(1);
00087 vFlux(0)+=data1.getMoles(1,0)/phasesVolume(1)*(1.0-mobW);
00088 vFlux(1)+=data1.getMoles(1,1)/phasesVolume(1)*(1.0-mobW);
00089 }
00090 else
00091 m_cl12=m_cl22=0;
00092
00093 m_flash.getPhasesVolume(P(face.cell2),data2,phasesVolume);
00094 m_Sw2=phasesVolume(0);
00095 mobW=m_mob(m_Sw2);
00096 m_mob.setParameters(m_vViscosities2);
00097
00098 if (phasesVolume(0) != 0.0)
00099 {
00100 m_cd11=data2.getMoles(0,0)/phasesVolume(0);
00101 m_cd21=data2.getMoles(0,1)/phasesVolume(0);
00102 vFlux(0)+=data2.getMoles(0,0)/phasesVolume(0)*mobW;
00103 vFlux(1)+=data2.getMoles(0,1)/phasesVolume(0)*mobW;
00104 }
00105 else
00106 m_cd11=m_cd21=0;
00107
00108 if (phasesVolume(1) != 0.0)
00109 {
00110 m_cd12=data2.getMoles(1,0)/phasesVolume(1);
00111 m_cd22=data2.getMoles(1,1)/phasesVolume(1);
00112 vFlux(0)+=data2.getMoles(1,0)/phasesVolume(1)*(1.0-mobW);
00113 vFlux(1)+=data2.getMoles(1,1)/phasesVolume(1)*(1.0-mobW);
00114 }
00115 else
00116 m_cd12=m_cd22=0;
00117
00118 vFlux*= m_vNormalVel(face.index)/2.0;
00119
00120 if (data1.getMoles(1,1) == 0.0 || data2.getMoles(1,1)==0.0)
00121 {
00122
00123 m_bUnsaturated=true;
00124 }
00125
00126
00127
00128 m_cell1=face.cell1;
00129 m_cell2=face.cell2;
00130 }
00131
00132
00133 double FaceFluxSimpleBlackOil::maxLocalCharVelocity(const FaceInfo &face,const VecDouble &Q1, const VecDouble &Q2)
00134 {
00135 assert(m_cell1==face.cell1);
00136 assert(m_cell2==face.cell2);
00137 double comp;
00138 if (m_bUnsaturated)
00139 comp=1.0;
00140 else
00141 comp=0.0;
00142
00143
00144 if (m_Sw1 > m_Sw2)
00145 std::swap(m_Sw1,m_Sw2);
00146
00147 if (m_vViscosities1(0)*m_vViscosities2(1) > m_vViscosities1(1)*m_vViscosities2(0))
00148 m_Dmob.setParameters(m_vViscosities1);
00149 else
00150 m_Dmob.setParameters(m_vViscosities2);
00151 m_Dmob.updateInflectionPoint();
00152
00153 return std::max(fabs(m_vNormalVel(face.index)*m_Dmob.getMaxNorm(m_Sw1,m_Sw2)),comp*m_vNormalVel(face.index));
00154
00155
00156 }
