henrylawobiblunt.cpp

Go to the documentation of this file.

#include "henrylawobiblunt.h"
#include "fluxforco2inj.h"
#include "numericmethods.h"
#include "arrayofvecdouble.h"

HenryLawObiBlunt::HenryLawObiBlunt(OrthoMesh &mesh,double Kh,double vn,double R,double T,double mc,double mv) 
  :FlashBase(2,2),
   m_Kh(Kh),
    m_c(-vn/(R*T)),
    m_mc(mc),
    m_mw(mv)
{
}




void HenryLawObiBlunt::execute() 
{
  TransportBase &trans=getTransportModule();
  //Get the necessary fields from the transport and dynamic module
  ArrayOfVecDouble &sol = trans.getSolutionAtCells();
  const VecDouble &vP = getDynamicModule().getPressureAtCells();
  
  assert(sol.size() == vP.size());


  //for each cell i in the mesh
  for (unsigned i=0;i<vP.size();i++)
  {
    double Cd = sol(i,FluxForCO2Inj::_Sd); //volume of dissolved CO2 per pore volume
    double Sc = sol(i,FluxForCO2Inj::_Sc); //Saturation of CO2
    double P = vP(i); //pression in the cell

    double Xd = (P/m_Kh)*exp(m_c*P); //Maximum molar fraction of dissolved co2
    if (Xd > 1.0)
    {
      printf("Houston we have a problem Xd == %g\n",Xd);
      throw new Exception("Finishing");
    }
    if (Sc >= 1.0)
      continue;
    
      
    //with Xd, we calculate the max amount of CO2 that can dissolve in water per pore volume
    double NDco2 = m_mw*(1.0-Sc - Cd)*Xd/(1.0-Xd);   


    double NTco2 = (Cd+Sc)*m_mc; //total moles of CO2 in the mixture
    double CdN,ScN;
    
    //Now compute the dissolution
    if (NTco2 >= NDco2)
    {
      //Ok, we have enough CO2 to dissolve into the water
      //dissolve it all.
      ScN = Sc + Cd - NDco2/m_mc; 
      CdN = NDco2/m_mc;
    }
    else
    {
      ScN = 0.0;
      CdN = (Cd +Sc);
    }
    sol(i,FluxForCO2Inj::_Sd)=CdN; //volume of dissolved CO2 per pore volume
    sol(i,FluxForCO2Inj::_Sc)=ScN; //Saturation of CO2
  } 
}
  
void HenryLawObiBlunt::updateDynamicModuleData()
{
}



 HenryLawObiBlunt::~HenryLawObiBlunt() 
{

}

---