flashsimpleblackoil.cpp

Go to the documentation of this file.

#include "flashsimpleblackoil.h"
#include "hdf5orthowriter.h"
#include "arrayofvecdouble.h"

//#define Rv 7.2518868898429352e-6
#define Rl 0.05
#define Cbl 2.31e-5

 FlashSimpleBlackOil::FlashSimpleBlackOil(OrthoMesh &mesh) 
   :FlashCompositional(2,2),m_mesh(mesh),m(2)
{
  compMass.reinit(2);
  compMass=1;
}



void FlashSimpleBlackOil::flash(double P, const VecDouble &compTotalMoles, FlashData &data) 
{
  assert(compTotalMoles.size()==2);
  assert(compTotalMoles(0)>=0.0);
  assert(compTotalMoles(1)>=0.0);
  assert(data.getNPhases()==2);
  assert(data.getNComponents()==2);
  
  data.allocateOwnMemory();

  //Find total moles.
  double no=compTotalMoles(OIL);
  double ng=compTotalMoles(GAS);

  //Find molar concentration of each component
  double ngo = Rl*P*no;
  

  data.setMoles(OIL,OIL,compTotalMoles(OIL));
  data.setMoles(GAS,OIL,0.0);
  

  if (ngo > ng)//Undersaturated case
  {
    //All gas is dissolved
    data.setMoles(OIL,GAS,ng);
    data.setMoles(GAS,GAS,0.0);
  }
  else
  {
    data.setMoles(OIL,GAS,ngo);
    data.setMoles(GAS,GAS,ng-ngo);
  }

}





 const VecDouble& FlashSimpleBlackOil::getComponentsMolarMass() const 
{
  return compMass;
}



 double FlashSimpleBlackOil::getFluidCompressibility(double P, FlashData &data) 
{
  double ngg = data.getMoles(GAS,GAS);
  if (ngg != 0.0) //Saturated case
  {
    double ng=data.getMoles(GAS,GAS)+data.getMoles(OIL,GAS);
    double no=data.getMoles(OIL,OIL);
    return -(-(0.3*no + 0.06*ng)/pow(6+0.06*P,2.0) + 0.0001*no);
  }
  else 
  {
    double no=data.getMoles(OIL,OIL);
    double pb=data.getMoles(OIL,GAS)/(Rl*no);
    return -(-(1+0.0001*pb)*Cbl*no/pow(1+Cbl*(P-pb),2));
  }


}

 void FlashSimpleBlackOil::getPhasesVolume(double P, const FlashData &data, VecDouble &phasesVol) 
{
  assert(data.getMoles(GAS,OIL)==0.0);
  assert(data.getNPhases()==2);
  assert(data.getNComponents()==2);
  
  //Sg = Bv*ngg == ngg/(6+0.06P)
  double ngg = data.getMoles(GAS,GAS);
  double noo = data.getMoles(OIL,OIL);
  
  if (ngg !=0.0) //Saturated case
  {
    phasesVol(GAS)=ngg/(6+0.06*P);
    phasesVol(OIL)=(1+1.e-4*P)*noo;
  }
  else
  {
    //No gas => no gas volume
    phasesVol(GAS)=0.0;

    //find Bubble pressure
    double pb = data.getMoles(OIL,GAS)/(Rl*noo);
    phasesVol(OIL)=(1+1.e-4*pb)*noo/(1+Cbl*(P-pb));
  }
}



 void FlashSimpleBlackOil::getPhasesViscosities(double P,const FlashData &data,VecDouble &visc) 
{
  assert(visc.size() == 2);
  visc(OIL)=1.10366e-12;
  visc(GAS)=9.19289e-14;
  return;
    


  visc(GAS)=0.012+3e-5*P;
  //Oil
  if (data.getMoles(GAS,GAS) != 0.0)//Saturated case
  {
    visc(OIL)=0.8-0.0001*P;
  }
  else
  {
    double pb= getBubblePressure(data);
    visc(OIL)=(0.8-0.0001*pb)*(1.0 + 0.0000678 * (P-pb));
  }

  //The formulas above ar for cp units. To change to PSI*day we need to multiply the viscosities by a factor
  visc*=1.6786775231481482e-12;
}


double FlashSimpleBlackOil::getBubblePressure(const FlashData &data) 
{
  assert(data.getMoles(GAS,OIL)==0.0);
  assert(data.getMoles(GAS,GAS)==0.0);
  return data.getMoles(OIL,GAS)/(Rl*data.getMoles(OIL,OIL));
}



 FlashSimpleBlackOil::~FlashSimpleBlackOil() 
{

}

void FlashSimpleBlackOil::flash(Index cell, FlashData &data) 
{

  getTransportModule().getSolutionAtCells().getVecValues(cell,&m);
  flash(getDynamicModule().getPressureAtCells()(cell),m,data);
}


void FlashSimpleBlackOil::getTotalVolumeDerivatives(double P,FlashData &data,VecDouble &dv_dm) 
{
  assert(dv_dm.size()  == 2);
  
  double ngg = data.getMoles(GAS,GAS);
  double no = data.getMoles(OIL,OIL);

  if (ngg != 0.0 )
  {
    //Saturated case
    
    dv_dm(OIL)=-0.05*P/(0.06*P + 6) + P/10000.0 + 1;
    dv_dm(GAS)=1.0/(0.06*P+6);
  }
  else //Unsaturated case
  {
    double ng=data.getMoles(OIL,GAS);
    dv_dm(OIL)=(((2.31e-5 *P + 1)*no +  -9.24e-4 * ng)*no -9.24e-7*ng*ng)/
      pow((2.31e-5*P+1)*no-ng*4.62e-4,2);
    dv_dm(GAS)=(4.62e-8*P + 0.002462)*no*no/pow( (2.31e-5*P + 1)*no - 4.62e-4*ng, 2);
  }

  
}


 void FlashSimpleBlackOil::printOutput() 
{
  HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter(); 
  unsigned nCells = m_mesh.numCells();
  VecDouble vc(nCells);
  VecDouble vcSatC(nCells);
  FlashData data(numPhases(),numComponents());
  VecDouble vv(m_mesh.numVertices());
  VecDouble phasesVol(2);
  for (unsigned i=0;i<nCells;i++)
  {
    flash(i,data);
    //vc(i)=data.getMoles(AQUEOUS,CO2);
    getPhasesVolume(getDynamicModule().getPressureAtCells()(i),data,phasesVol);
    vcSatC(i)=phasesVol(1)/(phasesVol(0)+phasesVol(1));
  }
  //m_mesh.projectCentralValuesAtVertices(vc,vv);
  //hdf5.writeScalarField(vv,"ngw");


  m_mesh.projectCentralValuesAtVertices(vcSatC,vv);
  hdf5.writeScalarField(vv,"SatC");

}

---