OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
MixtureBO3_ODGW.cpp
Go to the documentation of this file.
1 
12 #include "MixtureBO.hpp"
13 
15  // BOMixture_ODGW
17 
18 BOMixture_ODGW::BOMixture_ODGW(const ParamReservoir& rs_param, const USI& i)
19 {
20  BOMixtureInit(rs_param);
21 
22  PVTW.Setup(rs_param.PVTW_T.data[i]);
23  PVCO.Setup(rs_param.PVCO_T.data[i]);
24  PVDG.Setup(rs_param.PVDG_T.data[i]);
25 
26  data.resize(6, 0);
27  cdata.resize(6, 0);
28 }
29 
30 
31 void BOMixture_ODGW::InitFlash(const OCP_DBL& Pin, const OCP_DBL& Pbbin, const OCP_DBL& Tin,
32  const OCP_DBL* Sjin, const OCP_DBL& Vpore,
33  const OCP_DBL* Ziin)
34 {
35  for (USI j = 0; j < 3; j++) phaseExist[j] = false;
36  for (USI i = 0; i < 3; i++) Ni[i] = 0;
37  fill(xij.begin(), xij.end(), 0.0);
38 
39 
40  P = Pin;
41  S[1] = Sjin[1]; // gas saturation
42  S[2] = Sjin[2]; // water saturation
43 
44  // Water Property
45  PVTW.Eval_All(0, P, data, cdata);
46  OCP_DBL Pw0 = data[0];
47  OCP_DBL bw0 = data[1];
48  OCP_DBL cbw = data[2];
49  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
50  OCP_DBL bwp = -cbw * bw0;
51 
52  mu[2] = data[3];
53  rho[2] = std_RhoW / bw;
54  xi[2] = 1 / (CONV1 * bw);
55  Ni[2] = Vpore * S[2] * xi[2];
56 
57  USI phasecae;
58 
59  if (1 - S[1] - S[2] < TINY) {
60  if (S[1] < TINY)
61  phasecae = PHASE_W; // case 1 : water, no oil, no gas
62  else
63  phasecae = PHASE_GW; // case 2 : water, gas, no oil
64  } else if (S[1] < TINY)
65  phasecae = PHASE_OW; // case 3 : water, oil, no gas
66  else
67  phasecae = PHASE_ODGW; // case 4 : water, oil, gas
68 
69  switch (phasecae) {
70  case PHASE_W: {
71  // water
72  phaseExist[2] = true;
73  S[0] = 0;
74  S[1] = 0;
75  S[2] = 1;
76  xij[2 * 3 + 2] = 1;
77 
78  // hypothetical Oil property
79  PVCO.Eval_All(0, P, data, cdata);
80  OCP_DBL rs = data[1];
81  OCP_DBL bo = data[2];
82 
83  // hypothetical Gas property
84  PVDG.Eval_All(0, P, data, cdata);
85  OCP_DBL bg = data[1] * (CONV1 / 1000);
86 
87  // total
88  v[0] = 0;
89  v[1] = 0;
90  v[2] = CONV1 * Ni[2] * bw;
91  vf = v[2];
92  vfp = CONV1 * Ni[2] * bwp;
93  vfi[0] = CONV1 * bo - 1000 * rs * bg;
94  vfi[1] = 1000 * bg;
95  vfi[2] = CONV1 * bw;
96 
97  break;
98  }
99  case PHASE_GW: {
100  // water, gas
101  phaseExist[1] = true;
102  phaseExist[2] = true;
103  xij[1 * 3 + 1] = 1;
104  xij[2 * 3 + 2] = 1;
105 
106  // hypothetical Oil property
107  PVCO.Eval_All(0, P, data, cdata);
108  OCP_DBL rs = data[1];
109  OCP_DBL bo = data[2];
110 
111  // gas property
112  PVDG.Eval_All(0, P, data, cdata);
113  OCP_DBL bg = data[1] * (CONV1 / 1000);
114  OCP_DBL cbg = cdata[1] * (CONV1 / 1000);
115 
116  mu[1] = data[2];
117  xi[1] = 1 / 1000 / bg;
118  rho[1] = std_RhoG / bg;
119  Ni[1] = Vpore * S[1] * xi[1];
120 
121  v[0] = 0;
122  v[1] = 1000 * Ni[1] * bg; // Ni[0] = 0;
123  v[2] = CONV1 * Ni[2] * bw;
124  // total
125  vf = v[1] + v[2];
126  S[0] = 0;
127  S[1] = v[1] / vf;
128  S[2] = v[2] / vf;
129  vfp = 1000 * Ni[1] * cbg + CONV1 * Ni[2] * bwp;
130  vfi[0] = CONV1 * bo - 1000 * rs * bg;
131  vfi[1] = 1000 * bg;
132  vfi[2] = CONV1 * bw;
133 
134  break;
135  }
136  case PHASE_OW: {
137  // water, oil, unsaturated
138  phaseExist[0] = true;
139  phaseExist[2] = true;
140 
141  // oil property
142  OCP_DBL Pbb = Pbbin;
143  PVCO.Eval_All(0, Pbb, data, cdata);
144  OCP_DBL rs = data[1];
145  OCP_DBL bosat = data[2];
146  OCP_DBL muosat = data[3];
147  OCP_DBL cbosat = data[4];
148  OCP_DBL cmuosat = data[5];
149  OCP_DBL bo = bosat * (1 - cbosat * (P - Pbb));
150  OCP_DBL bop = -bosat * cbosat;
151  OCP_DBL dBo_drs = bo / bosat * cdata[2] +
152  bosat * (cdata[4] * (Pbb - P) + cbosat * cdata[0]);
153  dBo_drs /= cdata[1];
154 
155  Ni[0] = Vpore * (1 - S[1] - S[2]) / (CONV1 * bo);
156  Ni[1] = Ni[0] * rs;
157  xi[0] = (1 + rs) / (CONV1 * bo);
158  rho[0] = (std_RhoO + (1000 / CONV1) * rs * std_RhoG) / bo;
159  mu[0] = muosat * (1 + cmuosat * (P - Pbb));
160 
161  xij[0 * 3 + 0] = Ni[0] / (Ni[0] + Ni[1]);
162  xij[0 * 3 + 1] = 1 - xij[0 * 3 + 0];
163  xij[1 * 3 + 1] = 1;
164  xij[2 * 3 + 2] = 1;
165 
166  // total
167  v[0] = CONV1 * Ni[0] * bo;
168  v[2] = CONV1 * Ni[2] * bw;
169  vf = v[0] + v[2];
170  S[0] = v[0] / vf;
171  S[1] = 0;
172  S[2] = v[2] / vf;
173  vfp = CONV1 * (Ni[0] * bop + Ni[2] * bwp);
174  vfi[0] = CONV1 * (bo - dBo_drs * (Ni[1] / Ni[0]));
175  vfi[1] = CONV1 * dBo_drs;
176  vfi[2] = CONV1 * bw;
177 
178  break;
179  }
180  case PHASE_ODGW: {
181  phaseExist[0] = true;
182  phaseExist[1] = true;
183  phaseExist[2] = true;
184 
185  // oil property
186  PVCO.Eval_All(0, P, data, cdata);
187  OCP_DBL rs = data[1];
188  OCP_DBL bo = data[2];
189  OCP_DBL crs = cdata[1];
190  OCP_DBL cbosat = cdata[2];
191 
192  mu[0] = data[3];
193  Ni[0] = Vpore * (1 - S[1] - S[2]) / (CONV1 * bo);
194  xi[0] = (1 + rs) / (CONV1 * bo);
195  rho[0] = (std_RhoO + (1000 / CONV1) * rs * std_RhoG) / bo;
196 
197  // gas property
198  PVDG.Eval_All(0, P, data, cdata);
199  OCP_DBL bg = data[1] * (CONV1 / 1000);
200  OCP_DBL cbg = cdata[1] * (CONV1 / 1000);
201  Ni[1] = Vpore * S[1] / bg / 1000 + Ni[0] * rs;
202  xi[1] = 1 / CONV1 / data[1];
203  rho[1] = std_RhoG / bg;
204  mu[1] = data[2];
205 
206  xij[0 * 3 + 0] = 1 / (1 + rs);
207  xij[0 * 3 + 1] = 1 - xij[0 * 3 + 0];
208  xij[1 * 3 + 1] = 1;
209  xij[2 * 3 + 2] = 1;
210 
211  // total
212  v[0] = CONV1 * Ni[0] * bo;
213  v[1] = 1000 * (Ni[1] - rs * Ni[0]) * bg;
214  v[2] = CONV1 * Ni[2] * bw;
215 
216  vf = v[0] + v[1] + v[2];
217  S[0] = v[0] / vf;
218  S[1] = v[1] / vf;
219  S[2] = v[2] / vf;
220  vfp = CONV1 * Ni[0] * cbosat +
221  1000 * (-crs * Ni[0] * bg + (Ni[1] - rs * Ni[0]) * cbg) +
222  CONV1 * Ni[2] * bwp;
223  vfi[0] = CONV1 * bo - 1000 * rs * bg;
224  vfi[1] = 1000 * bg;
225  vfi[2] = CONV1 * bw;
226 
227  break;
228  }
229  }
230 }
231 
232 
233 
234 void BOMixture_ODGW::Flash(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Niin, const USI& ftype, const USI& lastNP,
235  const OCP_DBL* lastKs)
236 {
237  for (USI j = 0; j < 3; j++) phaseExist[j] = false;
238  fill(xij.begin(), xij.end(), 0.0);
239 
240  P = Pin;
241  Nt = 0;
242  for (USI i = 0; i < numCom; i++) {
243  Ni[i] = Niin[i];
244  Nt += Ni[i];
245  }
246 
247  // Water property
248  PVTW.Eval_All(0, P, data, cdata);
249  OCP_DBL Pw0 = data[0];
250  OCP_DBL bw0 = data[1];
251  OCP_DBL cbw = data[2];
252  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
253  OCP_DBL bwp = -cbw * bw0;
254 
255  mu[2] = data[3];
256  xi[2] = 1 / (CONV1 * bw);
257  rho[2] = std_RhoW / bw;
258 
259  USI phasecase;
260  OCP_DBL Rs_sat = PVCO.Eval(0, P, 1);
261 
262  if (Ni[0] < Nt * TINY) {
263  if (Ni[1] <= Ni[0] * Rs_sat)
264  phasecase = PHASE_W; // water, no oil, no gas
265  else
266  phasecase = PHASE_GW; // water, gas, no oil
267  } else if (Ni[1] <= Ni[0] * Rs_sat)
268  phasecase = PHASE_OW; // water, oil, no gas
269  else
270  phasecase = PHASE_ODGW; // water, oil ,gas
271 
272  switch (phasecase) {
273  case PHASE_W: {
274  // water
275  phaseExist[2] = true;
276  S[0] = 0;
277  S[1] = 0;
278  S[2] = 1;
279  xij[2 * 3 + 2] = 1;
280 
281  // hypothetical Oil property
282  PVCO.Eval_All(0, P, data, cdata);
283  OCP_DBL rs = data[1];
284  OCP_DBL bo = data[2];
285 
286  // hypothetical Gas property
287  PVDG.Eval_All(0, P, data, cdata);
288  OCP_DBL bg = data[1] * (CONV1 / 1000);
289 
290  // total
291  v[0] = 0;
292  v[1] = 0;
293  v[2] = CONV1 * Ni[2] * bw;
294  vf = v[2];
295  vfp = CONV1 * Ni[2] * bwp;
296  vfi[0] = CONV1 * bo - 1000 * rs * bg;
297  vfi[1] = 1000 * bg;
298  vfi[2] = CONV1 * bw;
299 
300  break;
301  }
302  case PHASE_GW: {
303  // water, gas
304  phaseExist[1] = true;
305  phaseExist[2] = true;
306  xij[1 * 3 + 1] = 1;
307  xij[2 * 3 + 2] = 1;
308 
309  // hypothetical Oil property
310  PVCO.Eval_All(0, P, data, cdata);
311  OCP_DBL rs = data[1];
312  OCP_DBL bo = data[2];
313 
314  // gas property
315  PVDG.Eval_All(0, P, data, cdata);
316  OCP_DBL bg = data[1] * (CONV1 / 1000);
317  OCP_DBL cbg = cdata[1] * (CONV1 / 1000);
318 
319  mu[1] = data[2];
320  xi[1] = 1 / 1000 / bg;
321  rho[1] = std_RhoG / bg;
322 
323  // total
324  v[0] = 0;
325  v[1] = 1000 * (Ni[1] - rs * Ni[0]) * bg;
326  v[2] = CONV1 * Ni[2] * bw;
327 
328 #ifdef _DEBUG
329  if (v[1] <= 0) {
330  OCP_ABORT("gas volume <= 0");
331  }
332 #endif // _DEBUG
333 
334  vf = v[1] + v[2];
335  S[0] = 0;
336  S[1] = v[1] / vf;
337  S[2] = v[2] / vf;
338  vfp = 1000 * Ni[1] * cbg + CONV1 * Ni[2] * bwp;
339  vfi[0] = CONV1 * bo - 1000 * rs * bg;
340  vfi[1] = 1000 * bg;
341  vfi[2] = CONV1 * bw;
342 
343  break;
344  }
345  case PHASE_OW: {
346  // water, oil
347  phaseExist[0] = true;
348  phaseExist[2] = true;
349  xij[0 * 3 + 0] = Ni[0] / (Ni[0] + Ni[1]);
350  xij[0 * 3 + 1] = 1 - xij[0 * 3 + 0];
351  xij[1 * 3 + 1] = 1;
352  xij[2 * 3 + 2] = 1;
353 
354  // oil property
355  OCP_DBL rs = Ni[1] / Ni[0];
356  PVCO.Eval_All(1, rs, data, cdata);
357  OCP_DBL pbb = data[0];
358  OCP_DBL bosat = data[2];
359  OCP_DBL muosat = data[3];
360  OCP_DBL cbosat = data[4];
361  OCP_DBL cmuosat = data[5];
362  OCP_DBL bo = bosat * (1 - cbosat * (P - pbb));
363  OCP_DBL bop = -bosat * cbosat;
364  OCP_DBL dBo_drs = bo / bosat * cdata[2] +
365  bosat * (cdata[4] * (pbb - P) + cbosat * cdata[0]);
366 
367  mu[0] = muosat * (1 + cmuosat * (P - pbb));
368  xi[0] = (1 + rs) / (CONV1 * bo);
369  rho[0] = (std_RhoO + (1000 / CONV1) * rs * std_RhoG) / bo;
370 
371  // total
372  v[0] = CONV1 * Ni[0] * bo;
373  v[2] = CONV1 * Ni[2] * bw;
374  vf = v[0] + v[2];
375  S[0] = v[0] / vf;
376  S[1] = 0;
377  S[2] = v[2] / vf;
378  vfp = CONV1 * (Ni[0] * bop + Ni[2] * bwp);
379  vfi[0] = CONV1 * (bo - dBo_drs * (Ni[1] / Ni[0]));
380  vfi[1] = CONV1 * dBo_drs;
381  vfi[2] = CONV1 * bw;
382 
383  break;
384  }
385  case PHASE_ODGW: {
386  phaseExist[0] = true;
387  phaseExist[1] = true;
388  phaseExist[2] = true;
389 
390  // oil property
391  PVCO.Eval_All(0, P, data, cdata);
392  OCP_DBL rs = data[1];
393  OCP_DBL bo = data[2];
394  OCP_DBL crs = cdata[1];
395  OCP_DBL cbosat = cdata[2];
396 
397  mu[0] = data[3];
398  xi[0] = (1 + rs) / (CONV1 * bo);
399  rho[0] = (std_RhoO + (1000 / CONV1) * rs * std_RhoG) / bo;
400 
401  // gas property
402  PVDG.Eval_All(0, P, data, cdata);
403  OCP_DBL bg = data[1] * (CONV1 / 1000);
404  OCP_DBL cbg = cdata[1] * (CONV1 / 1000);
405 
406  mu[1] = data[2];
407  xi[1] = 1 / CONV1 / data[1];
408  rho[1] = std_RhoG / bg;
409 
410  // total
411  xij[0 * 3 + 0] = 1 / (1 + rs);
412  xij[0 * 3 + 1] = 1 - xij[0 * 3 + 0];
413  xij[1 * 3 + 1] = 1;
414  xij[2 * 3 + 2] = 1;
415 
416  v[0] = CONV1 * Ni[0] * bo;
417  v[1] = 1000 * (Ni[1] - rs * Ni[0]) * bg;
418  v[2] = CONV1 * Ni[2] * bw;
419  vf = v[0] + v[1] + v[2];
420  S[0] = v[0] / vf;
421  S[1] = v[1] / vf;
422  S[2] = v[2] / vf;
423  vfp = CONV1 * Ni[0] * cbosat +
424  1000 * (-crs * Ni[0] * bg + (Ni[1] - rs * Ni[0]) * cbg) +
425  CONV1 * Ni[2] * bwp;
426  vfi[0] = CONV1 * bo - 1000 * rs * bg;
427  vfi[1] = 1000 * bg;
428  vfi[2] = CONV1 * bw;
429 
430  break;
431  }
432  }
433 }
434 
435 void BOMixture_ODGW::FlashDeriv(const OCP_DBL& Pin, const OCP_DBL& Tin,
436  const OCP_DBL* Niin, const USI& ftype, const USI& lastNP,
437  const OCP_DBL* lastKs)
438 {
439  for (USI j = 0; j < 3; j++) {
440  phaseExist[j] = false;
441  rhoP[j] = 0; xiP[j] = 0; muP[j] = 0;
442  }
443  fill(xij.begin(), xij.end(), 0.0);
444  fill(rhox.begin(), rhox.end(), 0.0);
445  fill(xix.begin(), xix.end(), 0.0);
446  fill(mux.begin(), mux.end(), 0.0);
447  fill(dXsdXp.begin(), dXsdXp.end(), 0.0);
448  fill(pEnumCom.begin(), pEnumCom.end(), 0.0);
449 
450  P = Pin;
451  Nt = 0;
452  for (USI i = 0; i < numCom; i++) {
453  Ni[i] = Niin[i];
454  Nt += Ni[i];
455  }
456 
457  // Water property
458  PVTW.Eval_All(0, P, data, cdata);
459  OCP_DBL Pw0 = data[0];
460  OCP_DBL bw0 = data[1];
461  OCP_DBL cbw = data[2];
462  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
463  OCP_DBL bwp = -cbw * bw0;
464 
465  mu[2] = data[3];
466  xi[2] = 1 / CONV1 / bw;
467  rho[2] = std_RhoW / bw;
468 
469  muP[2] = cdata[3];
470  xiP[2] = -bwp / (bw * bw * CONV1);
471  rhoP[2] = CONV1 * xiP[2] * std_RhoW;
472 
473  USI phasecase;
474  OCP_DBL Rs_sat = PVCO.Eval(0, P, 1);
475 
476  if (Ni[0] < Nt * TINY) {
477  if (Ni[1] <= Ni[0] * Rs_sat)
478  phasecase = PHASE_W; // water, no oil, no gas
479  else
480  phasecase = PHASE_GW; // water, gas, no oil
481  }
482  // Beacareful when NO = NG = 0, if it's possible.
483  else if (Ni[1] <= Ni[0] * Rs_sat)
484  phasecase = PHASE_OW; // water, oil, no gas
485  else
486  phasecase = PHASE_ODGW; // water, oil ,gas
487 
488  switch (phasecase) {
489  case PHASE_W: {
490  // water
491  phaseExist[2] = true;
492  S[0] = 0;
493  S[1] = 0;
494  S[2] = 1;
495  xij[2 * 3 + 2] = 1;
496 
497  // hypothetical Oil property
498  PVCO.Eval_All(0, P, data, cdata);
499  OCP_DBL rs = data[1];
500  OCP_DBL bo = data[2];
501 
502  // hypothetical Gas property
503  PVDG.Eval_All(0, P, data, cdata);
504  OCP_DBL bg = data[1] * (CONV1 / 1000);
505 
506  // total
507  v[0] = 0;
508  v[1] = 0;
509  v[2] = CONV1 * Ni[2] * bw;
510  vf = v[2];
511  vfp = CONV1 * Ni[2] * bwp;
512  vfi[0] = CONV1 * bo - 1000 * rs * bg;
513  vfi[1] = 1000 * bg;
514  vfi[2] = CONV1 * bw;
515 
516  // dXsdXp[0] = 0; // dSo / dP
517  dXsdXp[1] = CONV1 * bo / vf; // dSo / dNo
518  // dXsdXp[2] = 0; // dSo / dNg;
519  // dXsdXp[3] = 0; // dSo / dNw;
520 
521  // dXsdXp[1 * 4 + 0] = 0; // dSg / dP
522  // dXsdXp[1 * 4 + 1] = 0; // dSg / dNo
523  dXsdXp[1 * 4 + 2] = 1000 * bg / vf; // dSg / dNg
524  // dXsdXp[1 * 4 + 3] = 0; // dSg / dNw
525 
526  dXsdXp[2 * 4 + 0] = (Ni[2] * bwp * CONV1 - S[2] * vfp) / vf; // dSw / dP
527  dXsdXp[2 * 4 + 1] = -S[2] * vfi[0] / vf; // dSw / dNo
528  dXsdXp[2 * 4 + 2] = -S[2] * vfi[1] / vf; // dSw / dNg
529  dXsdXp[2 * 4 + 3] = (CONV1 * bw - S[2] * vfi[2]) / vf; // dSw / dNw
530 
531  //pEnumCom[0] = 0; pEnumCom[1] = 0; pEnumCom[2] = 0;
532 
533  break;
534  }
535  case PHASE_GW: {
536  // water, gas
537  phaseExist[1] = true;
538  phaseExist[2] = true;
539  xij[1 * 3 + 1] = 1;
540  xij[2 * 3 + 2] = 1;
541 
542  // hypothetical Oil property
543  PVCO.Eval_All(0, P, data, cdata);
544  OCP_DBL rs = data[1];
545  OCP_DBL bo = data[2];
546 
547  // gas property
548  PVDG.Eval_All(0, P, data, cdata);
549  OCP_DBL bg = data[1] * (CONV1 / 1000);
550  OCP_DBL cbg = cdata[1] * (CONV1 / 1000);
551 
552  mu[1] = data[2];
553  xi[1] = 1 / 1000 / bg;
554  rho[1] = std_RhoG / bg;
555 
556  muP[1] = cdata[2];
557  xiP[1] = -cdata[1] / (CONV1 * data[1] * data[1]);
558  rhoP[1] = 1000 * std_RhoG * xiP[1];
559 
560  // total
561  v[0] = 0;
562  v[1] = 1000 * (Ni[1] - rs * Ni[0]) * bg;
563  v[2] = CONV1 * Ni[2] * bw;
564 
565 #ifdef _DEBUG
566  if (v[1] <= 0) {
567  OCP_ABORT("gas volume <= 0");
568  }
569 #endif // _DEBUG
570 
571  vf = v[1] + v[2];
572  S[0] = 0;
573  S[1] = v[1] / vf;
574  S[2] = v[2] / vf;
575  vfp = 1000 * Ni[1] * cbg + CONV1 * Ni[2] * bwp;
576  vfi[0] = CONV1 * bo - 1000 * rs * bg;
577  vfi[1] = 1000 * bg;
578  vfi[2] = CONV1 * bw;
579 
580  dXsdXp[0] = 0; // dSo / dP
581  dXsdXp[1] = CONV1 * bo / vf; // dSo / dNo
582  dXsdXp[2] = 0; // dSo / dNg
583  dXsdXp[3] = 0; // dSo / dNw
584 
585  dXsdXp[1 * 4 + 0] = (1000 * Ni[1] * cbg - S[1] * vfp) / vf; // dSg / dP
586  dXsdXp[1 * 4 + 1] = (-1000 * rs * bg - S[1] * vfi[0]) / vf; // dSg / dNo
587  dXsdXp[1 * 4 + 2] = (1000 * bg - S[1] * vfi[1]) / vf; // dSg / dNg
588  dXsdXp[1 * 4 + 3] = -S[1] * vfi[2] / vf; // dSg / dNw
589 
590  dXsdXp[2 * 4 + 0] = (CONV1 * Ni[2] * bwp - S[2] * vfp) / vf; // dSw / dP
591  dXsdXp[2 * 4 + 1] = -S[2] * vfi[0] / vf; // dSw / dNo
592  dXsdXp[2 * 4 + 2] = -S[2] * vfi[1] / vf; // dSw / dNg
593  dXsdXp[2 * 4 + 3] = (CONV1 * bw - S[2] * vfi[2]) / vf; // dSw / dNw
594 
595  //pEnumCom[0] = 0; pEnumCom[1] = 0; pEnumCom[2] = 0;
596 
597  break;
598  }
599  case PHASE_OW: {
600  // water, oil
601  phaseExist[0] = true;
602  phaseExist[2] = true;
603  xij[0 * 3 + 0] = Ni[0] / (Ni[0] + Ni[1]);
604  xij[0 * 3 + 1] = 1 - xij[0 * 3 + 0];
605  xij[1 * 3 + 1] = 1;
606  xij[2 * 3 + 2] = 1;
607 
608  // oil property
609  OCP_DBL rs = Ni[1] / Ni[0];
610  PVCO.Eval_All(1, rs, data, cdata);
611  OCP_DBL pbb = data[0];
612  OCP_DBL bosat = data[2];
613  OCP_DBL muosat = data[3];
614  OCP_DBL cbosat = data[4];
615  OCP_DBL cmuosat = data[5];
616  OCP_DBL bo = bosat * (1 - cbosat * (P - pbb));
617  OCP_DBL bop = -bosat * cbosat;
618  OCP_DBL dBo_drs = (1 - cbosat * (P - pbb)) * cdata[2] +
619  bosat * (cdata[4] * (pbb - P) + cbosat * cdata[0]);
620 
621  mu[0] = muosat * (1 + cmuosat * (P - pbb));
622  xi[0] = (1 + rs) / (CONV1 * bo);
623  rho[0] = (std_RhoO + (1000 / CONV1) * rs * std_RhoG) / bo;
624 
625  muP[0] = muosat * cmuosat;
626  xiP[0] = -(1 + rs) * bop / (CONV1 * bo * bo);
627  rhoP[0] = -(std_RhoO + (1000/CONV1) * rs * std_RhoG) / (bo * bo) * bop;
628 
629  OCP_DBL muo_rs = mu[0] / muosat * cdata[3] + muosat * (cdata[5] * (P - pbb) - cmuosat * cdata[0]);
630  OCP_DBL xio_rs = 1 / CONV1 / bo - (1 + rs) * dBo_drs / (CONV1 * bo * bo);
631  OCP_DBL rhoo_rs = (1000 / CONV1) * std_RhoG / bo - (std_RhoO + (1000 / CONV1) * rs * std_RhoG) * dBo_drs / (bo * bo);
632 
633  // total
634  v[0] = CONV1 * Ni[0] * bo;
635  v[2] = CONV1 * Ni[2] * bw;
636  vf = v[0] + v[2];
637  S[0] = v[0] / vf;
638  S[1] = 0;
639  S[2] = v[2] / vf;
640  vfp = CONV1 * (Ni[0] * bop + Ni[2] * bwp);
641  vfi[0] = CONV1 * (bo - dBo_drs * (Ni[1] / Ni[0]));
642  vfi[1] = CONV1 * dBo_drs;
643  vfi[2] = CONV1 * bw;
644 
645 
646  dXsdXp[0] = (-CONV1 * Ni[0] * cbosat * bosat - S[0] * vfp) / vf; // dSo / dP
647  dXsdXp[1] = (CONV1 * bo - CONV1 * dBo_drs * (Ni[1] / Ni[0]) - S[0] * vfi[0]) / vf; // dSo / dNo
648  dXsdXp[2] = (CONV1 * dBo_drs - S[0] * vfi[1]) / vf; // dSo / dNg
649  dXsdXp[3] = -S[0] / vf * vfi[2]; // dSo / dNw
650 
651  dXsdXp[2 * 4 + 0] = (Ni[2] * bwp * CONV1 - S[2] * vfp) / vf; // dSw / dP
652  dXsdXp[2 * 4 + 1] = -S[2] * vfi[0] / vf; // dSw / dNo
653  dXsdXp[2 * 4 + 2] = -S[2] * vfi[1] / vf; // dSw / dNg
654  dXsdXp[2 * 4 + 3] = (CONV1 * bw - S[2] * vfi[2]) / vf; // dSw / dNw
655 
656  dXsdXp[3 * 4 + 1] = Ni[1] / pow((Ni[0] + Ni[1]), 2); // d Xoo / d No
657  dXsdXp[3 * 4 + 2] = -Ni[0] / pow((Ni[0] + Ni[1]), 2); // d Xoo / d Ng
658  dXsdXp[4 * 4 + 1] = -dXsdXp[3 * 4 + 1]; // d Xgo / d No
659  dXsdXp[4 * 4 + 2] = -dXsdXp[3 * 4 + 2]; // d Xgo / d Ng
660 
661  pEnumCom[0] = 2; pEnumCom[1] = 0; pEnumCom[2] = 0;
662 
663  OCP_DBL tmp = xij[0] * xij[0];
664 
665  mux[0] = -muo_rs * xij[1] / tmp; // dMuo / dXoo
666  mux[1] = muo_rs / xij[0]; // dMuo / dXgo
667 
668  xix[0] = -xio_rs * xij[1] / tmp; // dXio / dXoo
669  xix[1] = xio_rs / xij[0]; // dXio / dXgo
670 
671  rhox[0] = -rhoo_rs * xij[1] / tmp; // dRhoo / dXoo
672  rhox[1] = rhoo_rs / xij[0]; // dRhoo / dXgo
673 
674  // New! right I think
675  //OCP_DBL tmp_new = (1 + rs) * (1 + rs);
676 
677  //mux[0] = -muo_rs * tmp_new; // dMuo / dXoo
678  //mux[1] = muo_rs * tmp_new; // dMuo / dXgo
679 
680  //xix[0] = -xio_rs * tmp_new; // dXio / dXoo
681  //xix[1] = xio_rs * tmp_new; // dXio / dXgo
682 
683  //rhox[0] = -rhoo_rs * tmp_new; // dRhoo / dXoo
684  //rhox[1] = rhoo_rs * tmp_new; // dRhoo / dXgo
685 
686  break;
687  }
688  case PHASE_ODGW: {
689  phaseExist[0] = true;
690  phaseExist[1] = true;
691  phaseExist[2] = true;
692 
693  // oil property
694  PVCO.Eval_All(0, P, data, cdata);
695  OCP_DBL rs = data[1];
696  OCP_DBL bo = data[2];
697  OCP_DBL crs = cdata[1];
698  OCP_DBL cbosat = cdata[2];
699 
700  mu[0] = data[3];
701  xi[0] = (1 + rs) / (CONV1 * bo);
702  rho[0] = (std_RhoO + (1000 / CONV1) * rs * std_RhoG) / bo;
703 
704  muP[0] = cdata[3];
705  xiP[0] = (crs * bo - (1 + rs) * cbosat) / (bo * bo) / CONV1;
706  rhoP[0] = (1000 / CONV1) * std_RhoG * crs / bo - (std_RhoO + (1000/CONV1) * rs * std_RhoG) * cbosat / (bo * bo);
707 
708  // gas property
709  PVDG.Eval_All(0, P, data, cdata);
710  OCP_DBL bg = data[1] * (CONV1 / 1000);
711  OCP_DBL cbg = cdata[1] * (CONV1 / 1000);
712 
713  mu[1] = data[2];
714  xi[1] = 1 / CONV1 / data[1];
715  rho[1] = std_RhoG / bg;
716 
717  muP[1] = cdata[2];
718  xiP[1] = -cdata[1] / (CONV1 * data[1] * data[1]);
719  rhoP[1] = 1000 * std_RhoG * xiP[1];
720 
721  // total
722  xij[0 * 3 + 0] = 1 / (1 + rs);
723  xij[0 * 3 + 1] = 1 - xij[0 * 3 + 0];
724  xij[1 * 3 + 1] = 1;
725  xij[2 * 3 + 2] = 1;
726 
727  v[0] = CONV1 * Ni[0] * bo;
728  v[1] = 1000 * (Ni[1] - rs * Ni[0]) * bg;
729  v[2] = CONV1 * Ni[2] * bw;
730  vf = v[0] + v[1] + v[2];
731  S[0] = v[0] / vf;
732  S[1] = v[1] / vf;
733  S[2] = v[2] / vf;
734  vfp = CONV1 * Ni[0] * cbosat +
735  1000 * (-crs * Ni[0] * bg + (Ni[1] - rs * Ni[0]) * cbg) +
736  CONV1 * Ni[2] * bwp;
737  vfi[0] = CONV1 * bo - 1000 * rs * bg;
738  vfi[1] = 1000 * bg;
739  vfi[2] = CONV1 * bw;
740 
741 
742  dXsdXp[0] = (CONV1 * Ni[0] * cbosat - S[0] * vfp)/ vf; // dSo / dP
743  dXsdXp[1] = (CONV1 * bo - S[0] * vfi[0]) / vf; // dSo / dNo
744  dXsdXp[2] = -S[0] * vfi[1] / vf; // dSo / dNg
745  dXsdXp[3] = -S[0] * vfi[2] / vf; // dSo / dNw
746 
747  dXsdXp[1 * 4 + 0] = (1000 * (Ni[1] - rs * Ni[0]) * cbg - 1000 * Ni[0] * bg * crs - S[1] * vfp) / vf; // dSg / dP
748  dXsdXp[1 * 4 + 1] = (-1000 * rs * bg - S[1] * vfi[0]) / vf; // dSg / dNo
749  dXsdXp[1 * 4 + 2] = (1000 * bg - S[1] * vfi[1]) / vf; // dSg / dNg
750  dXsdXp[1 * 4 + 3] = -S[1] * vfi[2] / vf; // dSg / dNw
751 
752  dXsdXp[2 * 4 + 0] = (Ni[2] * bwp * CONV1 - S[2] * vfp) / vf; // dSw / dP
753  dXsdXp[2 * 4 + 1] = -S[2] * vfi[0] / vf; // dSw / dNo
754  dXsdXp[2 * 4 + 2] = -S[2] * vfi[1] / vf; // dSw / dNg
755  dXsdXp[2 * 4 + 3] = (CONV1 * bw - S[2] * vfi[2]) / vf; // dSw / dNw
756 
757  dXsdXp[3 * 4 + 0] = -crs / ((1 + rs) * (1 + rs)); // d Xoo / dP
758  dXsdXp[4 * 4 + 0] = -dXsdXp[3 * 4 + 0]; // d Xgo / dP
759 
760  pEnumCom[0] = 2; pEnumCom[1] = 0; pEnumCom[2] = 0;
761 
762  break;
763  }
764  }
765 
766 #ifdef Debug
767  for (auto v : dXsdXp) {
768  if (!isfinite(v)) {
769  OCP_ABORT("Nan in dXsdXp!");
770  }
771  }
772 #endif
773 }
774 
775 OCP_DBL BOMixture_ODGW::XiPhase(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Ziin)
776 {
777  if (Ziin[1] > 1 - TINY) {
778  // inj fluid is gas
779  OCP_DBL bg = PVDG.Eval(0, Pin, 1);
780  // OCP_DBL xig = 1 / (CONV1 * bg);
781  OCP_DBL xig = 1 / CONV1 / bg;
782  return xig;
783  } else if (Ziin[2] > 1 - TINY) {
784  // inj fluid is water
785 
786  PVTW.Eval_All(0, Pin, data, cdata);
787  OCP_DBL Pw0 = data[0];
788  OCP_DBL bw0 = data[1];
789  OCP_DBL cbw = data[2];
790  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
791  // OCP_DBL xiw = 1 / (CONV1 * bw);
792  OCP_DBL xiw = 1 / CONV1 / bw;
793  return xiw;
794  } else {
795  OCP_ABORT("Wrong Zi!");
796  }
797 }
798 
799 OCP_DBL BOMixture_ODGW::RhoPhase(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Ziin)
800 {
801  if (Ziin[1] > 1 - TINY) {
802  // inj fluid is gas
803  OCP_DBL bg = PVDG.Eval(0, Pin, 1);
804  OCP_DBL rhog = (1000 / CONV1) * std_RhoG / bg;
805  return rhog;
806  } else if (Ziin[2] > 1 - TINY) {
807  // inj fluid is water
808 
809  PVTW.Eval_All(0, Pin, data, cdata);
810  OCP_DBL Pw0 = data[0];
811  OCP_DBL bw0 = data[1];
812  OCP_DBL cbw = data[2];
813  OCP_DBL bw = bw0 * (1 - cbw * (P - Pw0));
814  OCP_DBL rhow = std_RhoW / bw;
815  return rhow;
816  } else {
817  OCP_ABORT("Wrong Zi!");
818  }
819 }
820 
821 OCP_DBL BOMixture_ODGW::GammaPhaseO(const OCP_DBL& Pin, const OCP_DBL& Pbbin)
822 {
823 
824  PVCO.Eval_All(0, Pbbin, data, cdata);
825  OCP_DBL rs = data[1];
826  OCP_DBL bosat = data[2];
827  OCP_DBL cbosat = data[4];
828  OCP_DBL bo = bosat * (1 - cbosat * (Pin - Pbbin));
829  OCP_DBL gammaO = (std_GammaO + (1000 / CONV1) * rs * std_GammaG) / bo;
830 
831  return gammaO;
832 }
833 
834 OCP_DBL BOMixture_ODGW::GammaPhaseG(const OCP_DBL& Pin)
835 {
836  if (PVDG.IsEmpty()) {
837  OCP_ABORT("PVDG is missing!");
838  }
839  OCP_DBL bg = (CONV1 / 1000) * PVDG.Eval(0, Pin, 1);
840 
841  return std_GammaG / bg;
842 }
843 
844 OCP_DBL BOMixture_ODGW::GammaPhaseW(const OCP_DBL& Pin)
845 {
846 
847  PVTW.Eval_All(0, Pin, data, cdata);
848  OCP_DBL Pw0 = data[0];
849  OCP_DBL bw0 = data[1];
850  OCP_DBL cbw = data[2];
851  OCP_DBL bw = (bw0 * (1 - cbw * (Pin - Pw0)));
852 
853  return std_GammaW / bw;
854 }
855 
856 
857 /*----------------------------------------------------------------------------*/
858 /* Brief Change History of This File */
859 /*----------------------------------------------------------------------------*/
860 /* Author Date Actions */
861 /*----------------------------------------------------------------------------*/
862 /* Shizhe Li Oct/01/2021 Create file */
863 /* Chensong Zhang Oct/15/2021 Format file */
864 /*----------------------------------------------------------------------------*/
MixtureBO.hpp
MixtureBO class declaration.
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:36
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
CONV1
const OCP_DBL CONV1
1 bbl = CONV1 ft3
Definition: OCPConst.hpp:59
PHASE_ODGW
const USI PHASE_ODGW
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:100
PHASE_W
const USI PHASE_W
Phase type = water only.
Definition: OCPConst.hpp:96
PHASE_OW
const USI PHASE_OW
Phase type = oil-water.
Definition: OCPConst.hpp:98
PHASE_GW
const USI PHASE_GW
Phase type = gas-water.
Definition: OCPConst.hpp:97
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
BOMixture_ODGW::InitFlash
void InitFlash(const OCP_DBL &Pin, const OCP_DBL &Pbbin, const OCP_DBL &Tin, const OCP_DBL *Sjin, const OCP_DBL &Vpore, const OCP_DBL *Ziin) override
flash calculation with saturation of phases.
Definition: MixtureBO3_ODGW.cpp:31
BOMixture_ODGW::GammaPhaseG
OCP_DBL GammaPhaseG(const OCP_DBL &Pin) override
return gamma of gas phase, gamma equals to mass density times gravity factor.
Definition: MixtureBO3_ODGW.cpp:834
BOMixture_ODGW::XiPhase
OCP_DBL XiPhase(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Ziin) override
Definition: MixtureBO3_ODGW.cpp:775
BOMixture_ODGW::Flash
void Flash(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const USI &ftype, const USI &lastNP, const OCP_DBL *lastKs) override
Flash calculation with moles of components.
Definition: MixtureBO3_ODGW.cpp:234
BOMixture_ODGW::FlashDeriv
void FlashDeriv(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Niin, const USI &ftype, const USI &lastNP, const OCP_DBL *lastKs) override
Flash calculation with moles of components and Calculate the derivative.
Definition: MixtureBO3_ODGW.cpp:435
BOMixture_ODGW::GammaPhaseW
OCP_DBL GammaPhaseW(const OCP_DBL &Pin) override
return gamma of water phase, gamma equals to mass density times gravity factor.
Definition: MixtureBO3_ODGW.cpp:844
BOMixture_ODGW::RhoPhase
OCP_DBL RhoPhase(const OCP_DBL &Pin, const OCP_DBL &Tin, const OCP_DBL *Ziin) override
return mass density of phase.
Definition: MixtureBO3_ODGW.cpp:799
BOMixture_ODGW::GammaPhaseO
OCP_DBL GammaPhaseO(const OCP_DBL &Pin, const OCP_DBL &Pbbin) override
return gamma of oil phase, gamma equals to mass density times gravity factor.
Definition: MixtureBO3_ODGW.cpp:821
BOMixture::std_GammaO
OCP_DBL std_GammaO
std_RhoO * gravity factor.
Definition: MixtureBO.hpp:76
BOMixture::std_RhoW
OCP_DBL std_RhoW
mass density of water phase in standard condition.
Definition: MixtureBO.hpp:79
BOMixture::std_RhoG
OCP_DBL std_RhoG
mass density of gas phase in standard condition.
Definition: MixtureBO.hpp:77
BOMixture::std_RhoO
OCP_DBL std_RhoO
< others.
Definition: MixtureBO.hpp:75
BOMixture::std_GammaW
OCP_DBL std_GammaW
std_RhoW * gravity factor.
Definition: MixtureBO.hpp:80
BOMixture::std_GammaG
OCP_DBL std_GammaG
std_RhoG * gravity factor.
Definition: MixtureBO.hpp:78
Mixture::rhoP
vector< OCP_DBL > rhoP
d rho / dP: numphase
Definition: Mixture.hpp:187
Mixture::pEnumCom
vector< USI > pEnumCom
see pEnumCom in bulk
Definition: Mixture.hpp:196
Mixture::numCom
USI numCom
num of components.
Definition: Mixture.hpp:157
Mixture::dXsdXp
vector< OCP_DBL > dXsdXp
the derivates of second variables wrt. primary variables
Definition: Mixture.hpp:195
Mixture::rho
vector< OCP_DBL > rho
mass density of phase: numPhase
Definition: Mixture.hpp:164
Mixture::mu
vector< OCP_DBL > mu
viscosity of phase: numPhase
Definition: Mixture.hpp:169
Mixture::phaseExist
vector< bool > phaseExist
existence of phase: numPhase
Definition: Mixture.hpp:162
Mixture::P
OCP_DBL P
pressure when flash calculation.
Definition: Mixture.hpp:158
Mixture::xix
vector< OCP_DBL > xix
d xi[j] / d x[i][j]: numphase * numCom
Definition: Mixture.hpp:192
Mixture::rhox
vector< OCP_DBL > rhox
d rho[j] / d x[i][j]: numphase * numCom
Definition: Mixture.hpp:193
Mixture::v
vector< OCP_DBL > v
volume of phase: numPhase;
Definition: Mixture.hpp:170
Mixture::S
vector< OCP_DBL > S
saturation of phase: numPhase
Definition: Mixture.hpp:163
Mixture::xij
vector< OCP_DBL > xij
Definition: Mixture.hpp:166
Mixture::Nt
OCP_DBL Nt
Total moles of Components.
Definition: Mixture.hpp:173
Mixture::Ni
vector< OCP_DBL > Ni
moles of component: numCom
Definition: Mixture.hpp:161
Mixture::vfi
vector< OCP_DBL > vfi
Definition: Mixture.hpp:182
Mixture::xiP
vector< OCP_DBL > xiP
d xi / dP: numphase
Definition: Mixture.hpp:186
Mixture::vf
OCP_DBL vf
volume of total fluids.
Definition: Mixture.hpp:172
Mixture::vfp
OCP_DBL vfp
Definition: Mixture.hpp:180
Mixture::muP
vector< OCP_DBL > muP
d mu / dP: numPhase
Definition: Mixture.hpp:185
Mixture::mux
vector< OCP_DBL > mux
d mu[j] / d x[i][j]: numphase * numCom
Definition: Mixture.hpp:191
Mixture::xi
vector< OCP_DBL > xi
molar density of phase: numPhase
Definition: Mixture.hpp:165
OCPTable::Eval_All
USI Eval_All(const USI &j, const OCP_DBL &val, vector< OCP_DBL > &outdata, vector< OCP_DBL > &slope)
Definition: OCPTable.cpp:53
OCPTable::Setup
void Setup(const vector< vector< OCP_DBL >> &src)
Setup tables from existing data of table.
Definition: OCPTable.cpp:32
OCPTable::IsEmpty
bool IsEmpty() const
judge if table is empty.
Definition: OCPTable.hpp:42
OCPTable::Eval
OCP_DBL Eval(const USI &j, const OCP_DBL &val, const USI &destj)
Definition: OCPTable.cpp:135
ParamReservoir::PVDG_T
TableSet PVDG_T
Table set of PVDG.
Definition: ParamReservoir.hpp:209
ParamReservoir::PVTW_T
TableSet PVTW_T
Table set of PVTW.
Definition: ParamReservoir.hpp:210
ParamReservoir::PVCO_T
TableSet PVCO_T
Table set of PVCO.
Definition: ParamReservoir.hpp:207
TableSet::data
vector< vector< vector< OCP_DBL > > > data
All table with the same name.
Definition: ParamReservoir.hpp:35
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1