23 BOMixtureInit(rs_param);
38 void BOMixture_ODGW::InitFlashIMPEC(
const OCP_DBL& Pin,
47 for (
USI i = 0; i < 3; i++)
Ni[i] = 0;
48 fill(
xij.begin(),
xij.end(), 0.0);
59 OCP_DBL bw = bw0 * (1 - cbw * (
P - Pw0));
65 Ni[2] = Vpore *
S[2] *
xi[2];
69 if (1 -
S[1] -
S[2] <
TINY) {
74 }
else if (
S[1] <
TINY)
104 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
129 xi[1] = 1 / 1000 / bg;
131 Ni[1] = Vpore *
S[1] *
xi[1];
134 vj[1] = 1000 *
Ni[1] * bg;
142 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
162 OCP_DBL bo = bosat * (1 - cbosat * (
P - Pbb));
164 OCP_DBL dBo_drs = bo / bosat * cdata[2] +
165 bosat * (cdata[4] * (Pbb -
P) + cbosat * cdata[0]);
168 Ni[0] = Vpore * (1 -
S[1] -
S[2]) / (
CONV1 * bo);
170 xi[0] = (1 + rs) / (
CONV1 * bo);
172 mu[0] = muosat * (1 + cmuosat * (
P - Pbb));
174 xij[0 * 3 + 0] =
Ni[0] / (
Ni[0] +
Ni[1]);
175 xij[0 * 3 + 1] = 1 -
xij[0 * 3 + 0];
207 Ni[0] = Vpore * (1 -
S[1] -
S[2]) / (
CONV1 * bo);
208 xi[0] = (1 + rs) / (
CONV1 * bo);
215 Ni[1] = Vpore *
S[1] / bg / 1000 +
Ni[0] * rs;
220 xij[0 * 3 + 0] = 1 / (1 + rs);
221 xij[0 * 3 + 1] = 1 -
xij[0 * 3 + 0];
227 vj[1] = 1000 * (
Ni[1] - rs *
Ni[0]) * bg;
235 1000 * (-crs *
Ni[0] * bg + (
Ni[1] - rs *
Ni[0]) * cbg) +
237 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
246 void BOMixture_ODGW::InitFlashFIM(
const OCP_DBL& Pin,
255 for (
USI i = 0; i < 3; i++)
Ni[i] = 0;
266 OCP_DBL bw = bw0 * (1 - cbw * (
P - Pw0));
269 Ni[2] = Vpore *
S[2] *
xi[2];
273 if (1 -
S[1] -
S[2] <
TINY) {
278 }
else if (
S[1] <
TINY)
296 xi[1] = 1 / 1000 / bg;
297 Ni[1] = Vpore *
S[1] *
xi[1];
309 OCP_DBL bo = bosat * (1 - cbosat * (
P - Pbb));
311 Ni[0] = Vpore * (1 -
S[1] -
S[2]) / (
CONV1 * bo);
322 Ni[0] = Vpore * (1 -
S[1] -
S[2]) / (
CONV1 * bo);
327 Ni[1] = Vpore *
S[1] / bg / 1000 +
Ni[0] * rs;
342 fill(
xij.begin(),
xij.end(), 0.0);
356 OCP_DBL bw = bw0 * (1 - cbw * (
P - Pw0));
367 if (
Ni[1] <=
Ni[0] * Rs_sat)
371 }
else if (
Ni[1] <=
Ni[0] * Rs_sat)
401 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
426 xi[1] = 1 / 1000 / bg;
431 vj[1] = 1000 * (
Ni[1] - rs *
Ni[0]) * bg;
443 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
454 xij[0 * 3 + 0] =
Ni[0] / (
Ni[0] +
Ni[1]);
455 xij[0 * 3 + 1] = 1 -
xij[0 * 3 + 0];
467 OCP_DBL bo = bosat * (1 - cbosat * (
P - pbb));
469 OCP_DBL dBo_drs = bo / bosat * cdata[2] +
470 bosat * (cdata[4] * (pbb -
P) + cbosat * cdata[0]);
472 mu[0] = muosat * (1 + cmuosat * (
P - pbb));
473 xi[0] = (1 + rs) / (
CONV1 * bo);
504 xi[0] = (1 + rs) / (
CONV1 * bo);
517 xij[0 * 3 + 0] = 1 / (1 + rs);
518 xij[0 * 3 + 1] = 1 -
xij[0 * 3 + 0];
523 vj[1] = 1000 * (
Ni[1] - rs *
Ni[0]) * bg;
530 1000 * (-crs *
Ni[0] * bg + (
Ni[1] - rs *
Ni[0]) * cbg) +
532 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
549 for (
USI j = 0; j < 3; j++) {
555 fill(
xij.begin(),
xij.end(), 0.0);
556 fill(
rhox.begin(),
rhox.end(), 0.0);
557 fill(
xix.begin(),
xix.end(), 0.0);
558 fill(
mux.begin(),
mux.end(), 0.0);
575 OCP_DBL bw = bw0 * (1 - cbw * (
P - Pw0));
590 if (
Ni[1] <=
Ni[0] * Rs_sat)
596 else if (
Ni[1] <=
Ni[0] * Rs_sat)
626 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
671 xi[1] = 1 / 1000 / bg;
675 xiP[1] = -cdata[1] / (
CONV1 * data[1] * data[1]);
680 vj[1] = 1000 * (
Ni[1] - rs *
Ni[0]) * bg;
692 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
702 dXsdXp[1 * 4 + 1] = (-1000 * rs * bg -
S[1] *
vfi[0]) /
vf;
723 xij[0 * 3 + 0] =
Ni[0] / (
Ni[0] +
Ni[1]);
724 xij[0 * 3 + 1] = 1 -
xij[0 * 3 + 0];
736 OCP_DBL bo = bosat * (1 - cbosat * (
P - pbb));
738 OCP_DBL dBo_drs = (1 - cbosat * (
P - pbb)) * cdata[2] +
739 bosat * (cdata[4] * (pbb -
P) + cbosat * cdata[0]);
741 mu[0] = muosat * (1 + cmuosat * (
P - pbb));
742 xi[0] = (1 + rs) / (
CONV1 * bo);
745 muP[0] = muosat * cmuosat;
746 xiP[0] = -(1 + rs) * bop / (
CONV1 * bo * bo);
750 OCP_DBL muo_rs =
mu[0] / muosat * cdata[3] +
751 muosat * (cdata[5] * (
P - pbb) - cmuosat * cdata[0]);
753 1 /
CONV1 / bo - (1 + rs) * dBo_drs / (
CONV1 * bo * bo);
807 mux[0] = -muo_rs *
xij[1] / tmp;
810 xix[0] = -xio_rs *
xij[1] / tmp;
813 rhox[0] = -rhoo_rs *
xij[1] / tmp;
844 xi[0] = (1 + rs) / (
CONV1 * bo);
848 xiP[0] = (crs * bo - (1 + rs) * cbosat) / (bo * bo) /
CONV1;
863 xiP[1] = -cdata[1] / (
CONV1 * data[1] * data[1]);
867 xij[0 * 3 + 0] = 1 / (1 + rs);
868 xij[0 * 3 + 1] = 1 -
xij[0 * 3 + 0];
873 vj[1] = 1000 * (
Ni[1] - rs *
Ni[0]) * bg;
880 1000 * (-crs *
Ni[0] * bg + (
Ni[1] - rs *
Ni[0]) * cbg) +
882 vfi[0] =
CONV1 * bo - 1000 * rs * bg;
891 dXsdXp[1 * 4 + 0] = (1000 * (
Ni[1] - rs *
Ni[0]) * cbg -
892 1000 *
Ni[0] * bg * crs -
S[1] *
vfP) /
894 dXsdXp[1 * 4 + 1] = (-1000 * rs * bg -
S[1] *
vfi[0]) /
vf;
903 dXsdXp[3 * 4 + 0] = -crs / ((1 + rs) * (1 + rs));
930 if (tarPhase ==
GAS) {
935 }
else if (tarPhase ==
WATER) {
941 OCP_DBL bw = bw0 * (1 - cbw * (Pin - Pw0));
957 if (tarPhase ==
OIL) {
958 PVCO.
Eval_All(0, Pbbin, data, cdata);
962 OCP_DBL bo = bosat * (1 - cbosat * (Pin - Pbbin));
965 }
else if (tarPhase ==
GAS) {
969 }
else if (tarPhase ==
WATER) {
974 OCP_DBL bw = bw0 * (1 - cbw * (Pin - Pw0));
983 const vector<SolventINJ>& sols,
987 const USI wellType = opt.WellType();
988 if (wellType ==
INJ) {
989 const string fluidName = opt.InjFluidType();
990 opt.SetInjFactor(1.0);
992 if (fluidName ==
"WAT") {
993 vector<OCP_DBL> tmpZi({0, 0, 1});
995 opt.SetInjProdPhase(
WATER);
996 }
else if (fluidName ==
"GAS") {
997 vector<OCP_DBL> tmpZi({0, 1, 0});
999 opt.SetInjProdPhase(
GAS);
1004 }
else if (wellType ==
PROD) {
1005 vector<OCP_DBL> tmpWght(3, 0);
1006 switch (opt.OptMode()) {
1019 tmpWght[0] = tmpWght[2] = 1;
1025 opt.SetProdPhaseWeight(tmpWght);
MixtureBO class declaration.
const OCP_DBL TINY
Small constant.
unsigned int USI
Generic unsigned integer.
const USI GAS
Fluid type = gas.
double OCP_DBL
Double precision.
const OCP_DBL CONV1
1 bbl = CONV1 ft3
const USI LRATE_MODE
Well option = fixed fluid rate???
const USI WRATE_MODE
Well option = fixed water rate.
const USI WATER
Fluid type = water.
const USI PHASE_ODGW
Phase type = oil-dry gas-water.
const USI ORATE_MODE
Well option = fixed oil rate.
const USI OIL
Fluid type = oil.
const USI BLKOIL_ODGW
black oil model with live oil, dry gas, water
const USI PHASE_W
Phase type = water only.
unsigned int OCP_USI
Long unsigned integer.
const USI PHASE_OW
Phase type = oil-water.
const USI PROD
Well type = producer.
const USI INJ
Well type = injector.
const USI BHP_MODE
Well option = fixed bottom-hole-pressure.
const USI GRATE_MODE
Well option = fixed gas rate.
const USI PHASE_GW
Phase type = gas-water.
#define OCP_ABORT(msg)
Abort if critical error happens.
void FlashFIM(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const OCP_DBL *Sjin, const USI &lastNP, const OCP_DBL *xijin, const OCP_USI &bId) override
Flash calculation with moles of components and Calculate the derivative.
OCP_DBL RhoPhase(const OCP_DBL &Pin, const OCP_DBL &Pbb, const OCP_DBL &Tin, const OCP_DBL *Ziin, const USI &tarPhase) override
return mass density of phase
void SetupWellOpt(WellOpt &opt, const vector< SolventINJ > &sols, const OCP_DBL &Psurf, const OCP_DBL &Tsurf) override
void FlashIMPEC(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const USI &lastNP, const OCP_DBL *xijin, const OCP_USI &bId) override
Flash calculation with moles of components.
OCP_DBL XiPhase(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Ziin, const USI &tarPhase) override
return mass density of phase
void Flash(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin) override
flash calculation with saturation of phases.
OCP_DBL std_RhoW
The density of water at surface conditions : lb/ft3.
OCP_DBL std_RhoG
The density of gas at surface conditions : lb/ft3.
OCP_DBL std_RhoO
< others.
vector< OCP_DBL > rhoP
d rho / dP: numphase
USI numCom
num of components.
vector< OCP_DBL > dXsdXp
derivatives of second variables wrt. primary variables
vector< OCP_DBL > rho
mass density of phase: numPhase
vector< OCP_DBL > mu
viscosity of phase: numPhase
OCP_DBL P
pressure when flash calculation.
vector< OCP_BOOL > pSderExist
Existence of derivative of phase saturation.
USI numPhase
num of phases.
vector< OCP_DBL > xix
d xi[j] / d x[i][j]: numphase * numCom
vector< OCP_DBL > vj
volume of phase: numPhase;
vector< OCP_DBL > rhox
d rho[j] / d x[i][j]: numphase * numCom
vector< OCP_DBL > S
saturation of phase: numPhase
vector< OCP_DBL > xij
Nij / nj: numPhase*numCom.
OCP_DBL Nt
Total moles of Components.
vector< OCP_DBL > Ni
moles of component: numCom
vector< OCP_BOOL > phaseExist
existence of phase: numPhase
vector< OCP_DBL > xiP
d xi / dP: numphase
OCP_DBL vf
volume of total fluids.
vector< USI > pVnumCom
num of variable components in the phase
vector< OCP_DBL > muP
d mu / dP: numPhase
vector< OCP_DBL > mux
d mu[j] / d x[i][j]: numphase * numCom
vector< OCP_DBL > xi
molar density of phase: numPhase
USI Eval_All(const USI &j, const OCP_DBL &val, vector< OCP_DBL > &outdata, vector< OCP_DBL > &slope)
Careful: the memory outdata and slope have not be allocated before.
void Setup(const vector< vector< OCP_DBL >> &src)
Setup tables from existing data of table.
OCP_DBL Eval(const USI &j, const OCP_DBL &val, const USI &destj)
TableSet PVDG_T
Table set of PVDG.
TableSet PVTW_T
Table set of PVTW.
TableSet PVCO_T
Table set of PVCO.
vector< vector< vector< OCP_DBL > > > data
All table with the same name.