My Project
Loading...
Searching...
No Matches
Brine_CO2.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_BINARY_COEFF_BRINE_CO2_HPP
29#define OPM_BINARY_COEFF_BRINE_CO2_HPP
30
33#include <opm/common/TimingMacros.hpp>
34
35#include <array>
36
37namespace Opm {
38namespace BinaryCoeff {
39
44template<class Scalar, class H2O, class CO2, bool verbose = true>
45class Brine_CO2 {
46 typedef ::Opm::IdealGas<Scalar> IdealGas;
47 static const int liquidPhaseIdx = 0; // index of the liquid phase
48 static const int gasPhaseIdx = 1; // index of the gas phase
49
50public:
58 template <class Evaluation>
59 static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure, bool extrapolate = false)
60 {
61 //Diffusion coefficient of water in the CO2 phase
62 Scalar k = 1.3806504e-23; // Boltzmann constant
63 Scalar c = 4; // slip parameter, can vary between 4 (slip condition) and 6 (stick condition)
64 Scalar R_h = 1.72e-10; // hydrodynamic radius of the solute
65 const Evaluation& mu = CO2::gasViscosity(temperature, pressure, extrapolate); // CO2 viscosity
66 return k / (c * M_PI * R_h) * (temperature / mu);
67 }
68
75 template <class Evaluation>
76 static Evaluation liquidDiffCoeff(const Evaluation& /*temperature*/, const Evaluation& /*pressure*/)
77 {
78 //Diffusion coefficient of CO2 in the brine phase
79 return 2e-9;
80 }
81
99 template <class Evaluation>
100 static void calculateMoleFractions(const Evaluation& temperature,
101 const Evaluation& pg,
102 const Evaluation& salinity,
103 const int knownPhaseIdx,
104 Evaluation& xlCO2,
105 Evaluation& ygH2O,
106 const int& activityModel,
107 bool extrapolate = false)
108 {
109 OPM_TIMEFUNCTION_LOCAL();
110
111 // Iterate or not?
112 bool iterate = false;
113 if ((activityModel == 1 && salinity > 0.0) || (activityModel == 2 && temperature > 372.15)) {
114 iterate = true;
115 }
116
117 // If both phases are present the mole fractions in each phase can be calculate with the mutual solubility
118 // function
119 if (knownPhaseIdx < 0) {
120 Evaluation molalityNaCl = massFracToMolality_(salinity); // mass fraction to molality of NaCl
121
122 // Duan-Sun model as given in Spycher & Pruess (2005) have a different fugacity coefficient formula and
123 // activity coefficient definition (not a true activity coefficient but a ratio).
124 // Technically only valid below T = 100 C, but we use low-temp. parameters and formulas even above 100 C as
125 // an approximation.
126 if (activityModel == 3) {
127 auto [xCO2, yH2O] = mutualSolubilitySpycherPruess2005_(temperature, pg, molalityNaCl, extrapolate);
128 xlCO2 = xCO2;
129 ygH2O = yH2O;
130
131 }
132 else {
133 // Fixed-point iterations to calculate solubility
134 if (iterate) {
135 auto [xCO2, yH2O] = fixPointIterSolubility_(temperature, pg, molalityNaCl, activityModel, extrapolate);
136 xlCO2 = xCO2;
137 ygH2O = yH2O;
138 }
139
140 // Solve mutual solubility equation with back substitution (no need for iterations)
141 else {
142 auto [xCO2, yH2O] = nonIterSolubility_(temperature, pg, molalityNaCl, activityModel, extrapolate);
143 xlCO2 = xCO2;
144 ygH2O = yH2O;
145 }
146 }
147 }
148
149 // if only liquid phase is present the mole fraction of CO2 in brine is given and
150 // and the virtual equilibrium mole fraction of water in the non-existing gas phase can be estimated
151 // with the mutual solubility function
152 else if (knownPhaseIdx == liquidPhaseIdx && activityModel == 3) {
153 Evaluation x_NaCl = salinityToMolFrac_(salinity);
154 const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
155 ygH2O = A * (1 - xlCO2 - x_NaCl);
156 }
157
158 // if only gas phase is present the mole fraction of water in the gas phase is given and
159 // and the virtual equilibrium mole fraction of CO2 in the non-existing liquid phase can be estimated
160 // with the mutual solubility function
161 else if (knownPhaseIdx == gasPhaseIdx && activityModel == 3) {
162 //y_H2o = fluidstate.
163 Evaluation x_NaCl = salinityToMolFrac_(salinity);
164 const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
165 xlCO2 = 1 - x_NaCl - ygH2O / A;
166 }
167 }
168
172 template <class Evaluation>
173 static Evaluation henry(const Evaluation& temperature, bool extrapolate = false)
174 { return fugacityCoefficientCO2(temperature, /*pressure=*/1e5, extrapolate)*1e5; }
175
184 template <class Evaluation>
185 static Evaluation fugacityCoefficientCO2(const Evaluation& temperature,
186 const Evaluation& pg,
187 const Evaluation& yH2O,
188 const bool highTemp,
189 bool extrapolate = false,
190 bool spycherPruess2005 = false)
191 {
192 OPM_TIMEFUNCTION_LOCAL();
193 Valgrind::CheckDefined(temperature);
194 Valgrind::CheckDefined(pg);
195
196 Evaluation V = 1 / (CO2::gasDensity(temperature, pg, extrapolate) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
197 Evaluation pg_bar = pg / 1.e5; // gas phase pressure in bar
198 Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
199
200 // Parameters in Redlich-Kwong equation
201 Evaluation a_CO2 = aCO2_(temperature, highTemp);
202 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
203 Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
204 Scalar b_CO2 = bCO2_(highTemp);
205 Evaluation b_mix = bMix_(yH2O, highTemp);
206
207 Evaluation lnPhiCO2;
208 if (spycherPruess2005) {
209 lnPhiCO2 = log(V / (V - b_CO2));
210 lnPhiCO2 += b_CO2 / (V - b_CO2);
211 lnPhiCO2 -= 2 * a_CO2 / (R * pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
212 lnPhiCO2 +=
213 a_CO2 * b_CO2
214 / (R
215 * pow(temperature, 1.5)
216 * b_CO2
217 * b_CO2)
218 * (log((V + b_CO2) / V)
219 - b_CO2 / (V + b_CO2));
220 lnPhiCO2 -= log(pg_bar * V / (R * temperature));
221 }
222 else {
223 lnPhiCO2 = (b_CO2 / b_mix) * (pg_bar * V / (R * temperature) - 1);
224 lnPhiCO2 -= log(pg_bar * (V - b_mix) / (R * temperature));
225 lnPhiCO2 += (2 * (yH2O * a_CO2_H2O + (1 - yH2O) * a_CO2) / a_mix - (b_CO2 / b_mix)) *
226 a_mix / (b_mix * R * pow(temperature, 1.5)) * log(V / (V + b_mix));
227 }
228 return exp(lnPhiCO2); // fugacity coefficient of CO2
229 }
230
239 template <class Evaluation>
240 static Evaluation fugacityCoefficientH2O(const Evaluation& temperature,
241 const Evaluation& pg,
242 const Evaluation& yH2O,
243 const bool highTemp,
244 bool extrapolate = false,
245 bool spycherPruess2005 = false)
246 {
247 OPM_TIMEFUNCTION_LOCAL();
248 Valgrind::CheckDefined(temperature);
249 Valgrind::CheckDefined(pg);
250
251 const Evaluation& V = 1 / (CO2::gasDensity(temperature, pg, extrapolate) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
252 const Evaluation& pg_bar = pg / 1.e5; // gas phase pressure in bar
253 Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
254
255 // Mixture parameter of Redlich-Kwong equation
256 Evaluation a_H2O = aH2O_(temperature, highTemp);
257 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
258 Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
259 Scalar b_H2O = bH2O_(highTemp);
260 Evaluation b_mix = bMix_(yH2O, highTemp);
261
262 Evaluation lnPhiH2O;
263 if (spycherPruess2005) {
264 lnPhiH2O =
265 log(V/(V - b_mix))
266 + b_H2O/(V - b_mix) - 2*a_CO2_H2O
267 / (R*pow(temperature, 1.5)*b_mix)*log((V + b_mix)/V)
268 + a_mix*b_H2O/(R*pow(temperature, 1.5)*b_mix*b_mix)
269 *(log((V + b_mix)/V) - b_mix/(V + b_mix))
270 - log(pg_bar*V/(R*temperature));
271 }
272 else {
273 lnPhiH2O = (b_H2O / b_mix) * (pg_bar * V / (R * temperature) - 1);
274 lnPhiH2O -= log(pg_bar * (V - b_mix) / (R * temperature));
275 lnPhiH2O += (2 * (yH2O * a_H2O + (1 - yH2O) * a_CO2_H2O) / a_mix - (b_H2O / b_mix)) *
276 a_mix / (b_mix * R * pow(temperature, 1.5)) * log(V / (V + b_mix));
277 }
278 return exp(lnPhiH2O); // fugacity coefficient of H2O
279 }
280
281private:
285 template <class Evaluation>
286 static Evaluation aCO2_(const Evaluation& temperature, const bool& highTemp)
287 {
288 if (highTemp) {
289 return 8.008e7 - 4.984e4 * temperature;
290 }
291 else {
292 return 7.54e7 - 4.13e4 * temperature;
293 }
294 }
295
299 template <class Evaluation>
300 static Evaluation aH2O_(const Evaluation& temperature, const bool& highTemp)
301 {
302 if (highTemp) {
303 return 1.337e8 - 1.4e4 * temperature;
304 }
305 else {
306 return 0.0;
307 }
308 }
309
313 template <class Evaluation>
314 static Evaluation aCO2_H2O_(const Evaluation& temperature, const Evaluation& yH2O, const bool& highTemp)
315 {
316 if (highTemp) {
317 // Pure parameters
318 Evaluation aCO2 = aCO2_(temperature, highTemp);
319 Evaluation aH2O = aH2O_(temperature, highTemp);
320
321 // Mixture Eq. (A-6)
322 Evaluation K_CO2_H2O = 0.4228 - 7.422e-4 * temperature;
323 Evaluation K_H2O_CO2 = 1.427e-2 - 4.037e-4 * temperature;
324 Evaluation k_CO2_H2O = yH2O * K_H2O_CO2 + (1 - yH2O) * K_CO2_H2O;
325
326 // Eq. (A-5)
327 return sqrt(aCO2 * aH2O) * (1 - k_CO2_H2O);
328 }
329 else {
330 return 7.89e7;
331 }
332 }
333
337 template <class Evaluation>
338 static Evaluation aMix_(const Evaluation& temperature, const Evaluation& yH2O, const bool& highTemp)
339 {
340 if (highTemp) {
341 // Parameters
342 Evaluation aCO2 = aCO2_(temperature, highTemp);
343 Evaluation aH2O = aH2O_(temperature, highTemp);
344 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
345
346 return yH2O * yH2O * aH2O + 2 * yH2O * (1 - yH2O) * a_CO2_H2O + (1 - yH2O) * (1 - yH2O) * aCO2;
347 }
348 else {
349 return aCO2_(temperature, highTemp);
350 }
351 }
352
356 static Scalar bCO2_(const bool& highTemp)
357 {
358 if (highTemp) {
359 return 28.25;
360 }
361 else {
362 return 27.8;
363 }
364 }
365
369 static Scalar bH2O_(const bool& highTemp)
370 {
371 if (highTemp) {
372 return 15.7;
373 }
374 else {
375 return 18.18;
376 }
377 }
378
382 template <class Evaluation>
383 static Evaluation bMix_(const Evaluation& yH2O, const bool& highTemp)
384 {
385 if (highTemp) {
386 // Parameters
387 Scalar bCO2 = bCO2_(highTemp);
388 Scalar bH2O = bH2O_(highTemp);
389
390 return yH2O * bH2O + (1 - yH2O) * bCO2;
391 }
392 else {
393 return bCO2_(highTemp);
394 }
395 }
396
400 template <class Evaluation>
401 static Evaluation V_avg_CO2_(const Evaluation& temperature, const bool& highTemp)
402 {
403 if (highTemp && (temperature > 373.15)) {
404 return 32.6 + 3.413e-2 * (temperature - 373.15);
405 }
406 else {
407 return 32.6;
408 }
409 }
410
414 template <class Evaluation>
415 static Evaluation V_avg_H2O_(const Evaluation& temperature, const bool& highTemp)
416 {
417 if (highTemp && (temperature > 373.15)) {
418 return 18.1 + 3.137e-2 * (temperature - 373.15);
419 }
420 else {
421 return 18.1;
422 }
423 }
424
428 template <class Evaluation>
429 static Evaluation AM_(const Evaluation& temperature, const bool& highTemp)
430 {
431 if (highTemp && temperature > 373.15) {
432 Evaluation deltaTk = temperature - 373.15;
433 return deltaTk * (-3.084e-2 + 1.927e-5 * deltaTk);
434 }
435 else {
436 return 0.0;
437 }
438 }
439
443 template <class Evaluation>
444 static Evaluation Pref_(const Evaluation& temperature, const bool& highTemp)
445 {
446 if (highTemp && temperature > 373.15) {
447 const Evaluation& temperatureCelcius = temperature - 273.15;
448 static const Scalar c[5] = { -1.9906e-1, 2.0471e-3, 1.0152e-4, -1.4234e-6, 1.4168e-8 };
449 return c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
450 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
451 }
452 else {
453 return 1.0;
454 }
455 }
456
460 template <class Evaluation>
461 static Evaluation activityCoefficientCO2_(const Evaluation& temperature,
462 const Evaluation& xCO2,
463 const bool& highTemp)
464 {
465 if (highTemp) {
466 // Eq. (13)
467 Evaluation AM = AM_(temperature, highTemp);
468 Evaluation lnGammaCO2 = 2 * AM * xCO2 * (1 - xCO2) * (1 - xCO2);
469 return exp(lnGammaCO2);
470 }
471 else {
472 return 1.0;
473 }
474 }
475
479 template <class Evaluation>
480 static Evaluation activityCoefficientH2O_(const Evaluation& temperature,
481 const Evaluation& xCO2,
482 const bool& highTemp)
483 {
484 if (highTemp) {
485 // Eq. (12)
486 Evaluation AM = AM_(temperature, highTemp);
487 Evaluation lnGammaH2O = (1 - 2 * (1 - xCO2)) * AM * xCO2 * xCO2;
488 return exp(lnGammaH2O);
489 }
490 else {
491 return 1.0;
492 }
493 }
494
500 template <class Evaluation>
501 static Evaluation salinityToMolFrac_(const Evaluation& salinity) {
502 OPM_TIMEFUNCTION_LOCAL();
503 const Scalar Mw = H2O::molarMass(); /* molecular weight of water [kg/mol] */
504 const Scalar Ms = 58.44e-3; /* molecular weight of NaCl [kg/mol] */
505
506 const Evaluation X_NaCl = salinity;
507 /* salinity: conversion from mass fraction to mol fraction */
508 const Evaluation x_NaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
509 return x_NaCl;
510 }
511
517#if 0
518 template <class Evaluation>
519 static Evaluation moleFracToMolality_(const Evaluation& x_NaCl)
520 {
521 // conversion from mol fraction to molality (dissolved CO2 neglected)
522 return 55.508 * x_NaCl / (1 - x_NaCl);
523 }
524#endif
525
526 template <class Evaluation>
527 static Evaluation massFracToMolality_(const Evaluation& X_NaCl)
528 {
529 const Scalar MmNaCl = 58.44e-3;
530 return X_NaCl / (MmNaCl * (1 - X_NaCl));
531 }
532
538 template <class Evaluation>
539 static Evaluation molalityToMoleFrac_(const Evaluation& m_NaCl)
540 {
541 // conversion from molality to mole fractio (dissolved CO2 neglected)
542 return m_NaCl / (55.508 + m_NaCl);
543 }
544
548 template <class Evaluation>
549 static std::pair<Evaluation, Evaluation> fixPointIterSolubility_(const Evaluation& temperature,
550 const Evaluation& pg,
551 const Evaluation& m_NaCl,
552 const int& activityModel,
553 bool extrapolate = false)
554 {
555 OPM_TIMEFUNCTION_LOCAL();
556 // Start point for fixed-point iterations as recommended below in section 2.2
557 Evaluation yH2O = H2O::vaporPressure(temperature) / pg; // ideal mixing
558 Evaluation xCO2 = 0.009; // same as ~0.5 mol/kg
559 Evaluation gammaNaCl = 1.0; // default salt activity coeff = 1.0
560
561 // We can pre-calculate Duan-Sun, Spycher & Pruess (2009) salt activity coeff.
562 if (m_NaCl > 0.0 && activityModel == 2) {
563 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
564 }
565
566 // Options
567 int max_iter = 100;
568 Scalar tol = 1e-8;
569 bool highTemp = true;
570 if (activityModel == 1) {
571 highTemp = false;
572 }
573 const bool iterate = true;
574
575 // Fixed-point loop x_i+1 = F(x_i)
576 for (int i = 0; i < max_iter; ++i) {
577 // Calculate activity coefficient for Rumpf et al (1994) model
578 if (m_NaCl > 0.0 && activityModel == 1) {
579 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, xCO2, activityModel);
580 }
581
582 // F(x_i) is the mutual solubilities
583 auto [xCO2_new, yH2O_new] = mutualSolubility_(temperature, pg, xCO2, yH2O, m_NaCl, gammaNaCl, highTemp,
584 iterate, extrapolate);
585
586 // Check for convergence
587 if (abs(xCO2_new - xCO2) < tol && abs(yH2O_new - yH2O) < tol) {
588 xCO2 = xCO2_new;
589 yH2O = yH2O_new;
590 break;
591 }
592
593 // Else update mole fractions for next iteration
594 else {
595 xCO2 = xCO2_new;
596 yH2O = yH2O_new;
597 }
598 }
599
600 return {xCO2, yH2O};
601 }
602
606 template <class Evaluation>
607 static std::pair<Evaluation, Evaluation> nonIterSolubility_(const Evaluation& temperature,
608 const Evaluation& pg,
609 const Evaluation& m_NaCl,
610 const int& activityModel,
611 bool extrapolate = false)
612 {
613 // Calculate activity coefficient for salt
614 Evaluation gammaNaCl = 1.0;
615 if (m_NaCl > 0.0 && activityModel > 0 && activityModel < 3) {
616 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
617 }
618
619 // Calculate mutual solubility.
620 // Note that we don't use xCO2 and yH2O input in low-temperature case, so we set them to 0.0
621 const bool highTemp = false;
622 const bool iterate = false;
623 auto [xCO2, yH2O] = mutualSolubility_(temperature, pg, Evaluation(0.0), Evaluation(0.0), m_NaCl, gammaNaCl,
624 highTemp, iterate, extrapolate);
625
626 return {xCO2, yH2O};
627 }
628
632 template <class Evaluation>
633 static std::pair<Evaluation, Evaluation> mutualSolubility_(const Evaluation& temperature,
634 const Evaluation& pg,
635 const Evaluation& xCO2,
636 const Evaluation& yH2O,
637 const Evaluation& m_NaCl,
638 const Evaluation& gammaNaCl,
639 const bool& highTemp,
640 const bool& iterate,
641 bool extrapolate = false)
642 {
643 // Calculate A and B (without salt effect); Eqs. (8) and (9)
644 const Evaluation& A = computeA_(temperature, pg, yH2O, xCO2, highTemp, extrapolate);
645 Evaluation B = computeB_(temperature, pg, yH2O, xCO2, highTemp, extrapolate);
646
647 // Add salt effect to B, Eq. (17)
648 B /= gammaNaCl;
649
650 // Compute yH2O and xCO2, Eqs. (B-7) and (B-2)
651 Evaluation yH2O_new = (1. - B) * 55.508 / ((1. / A - B) * (2 * m_NaCl + 55.508) + 2 * m_NaCl * B);
652 Evaluation xCO2_new;
653 if (iterate) {
654 xCO2_new = B * (1 - yH2O);
655 }
656 else {
657 xCO2_new = B * (1 - yH2O_new);
658 }
659
660 return {xCO2_new, yH2O_new};
661 }
662
666 template <class Evaluation>
667 static std::pair<Evaluation, Evaluation> mutualSolubilitySpycherPruess2005_(const Evaluation& temperature,
668 const Evaluation& pg,
669 const Evaluation& m_NaCl,
670 bool extrapolate = false)
671 {
672 // Calculate A and B (without salt effect); Eqs. (8) and (9)
673 const Evaluation& A = computeA_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
674 const Evaluation& B = computeB_(temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
675
676 // Mole fractions and molality in pure water
677 Evaluation yH2O = (1 - B) / (1. / A - B);
678 Evaluation xCO2 = B * (1 - yH2O);
679
680 // Modifiy mole fractions with Duan-Sun "activity coefficient" if salt is involved
681 if (m_NaCl > 0.0) {
682 const Evaluation& gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), 3);
683
684 // Molality with salt
685 Evaluation mCO2 = (xCO2 * 55.508) / (1 - xCO2); // pure water
686 mCO2 /= gammaNaCl;
687 xCO2 = mCO2 / (m_NaCl + 55.508 + mCO2);
688
689 // new yH2O with salt
690 const Evaluation& xNaCl = molalityToMoleFrac_(m_NaCl);
691 yH2O = A * (1 - xCO2 - xNaCl);
692 }
693
694 return {xCO2, yH2O};
695 }
696
705 template <class Evaluation>
706 static Evaluation computeA_(const Evaluation& temperature,
707 const Evaluation& pg,
708 const Evaluation& yH2O,
709 const Evaluation& xCO2,
710 const bool& highTemp,
711 bool extrapolate = false,
712 bool spycherPruess2005 = false)
713 {
714 OPM_TIMEFUNCTION_LOCAL();
715 // Intermediate calculations
716 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp); // pressure range [bar] from pref to pg[bar]
717 Evaluation v_av_H2O = V_avg_H2O_(temperature, highTemp); // average partial molar volume of H2O [cm^3/mol]
718 Evaluation k0_H2O = equilibriumConstantH2O_(temperature, highTemp); // equilibrium constant for H2O at 1 bar
719 Evaluation phi_H2O = fugacityCoefficientH2O(temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005); // fugacity coefficient of H2O for the water-CO2 system
720 Evaluation gammaH2O = activityCoefficientH2O_(temperature, xCO2, highTemp);
721
722 // In the intermediate temperature range 99-109 C, equilibrium constants and fugacity coeff. are linearly
723 // weighted
724 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
725 const Evaluation weight = (382.15 - temperature) / 10.;
726 const Evaluation& k0_H2O_low = equilibriumConstantH2O_(temperature, false);
727 const Evaluation& phi_H2O_low = fugacityCoefficientH2O(temperature, pg, Evaluation(0.0), false, extrapolate);
728 k0_H2O = k0_H2O * (1 - weight) + k0_H2O_low * weight;
729 phi_H2O = phi_H2O * (1 - weight) + phi_H2O_low * weight;
730 }
731
732 // Eq. (10)
733 const Evaluation& pg_bar = pg / 1.e5;
734 Scalar R = IdealGas::R * 10;
735 return k0_H2O * gammaH2O / (phi_H2O * pg_bar) * exp(deltaP * v_av_H2O / (R * temperature));
736 }
737
746 template <class Evaluation>
747 static Evaluation computeB_(const Evaluation& temperature,
748 const Evaluation& pg,
749 const Evaluation& yH2O,
750 const Evaluation& xCO2,
751 const bool& highTemp,
752 bool extrapolate = false,
753 bool spycherPruess2005 = false)
754 {
755 OPM_TIMEFUNCTION_LOCAL();
756 // Intermediate calculations
757 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp); // pressure range [bar] from pref to pg[bar]
758 Evaluation v_av_CO2 = V_avg_CO2_(temperature, highTemp); // average partial molar volume of CO2 [cm^3/mol]
759 Evaluation k0_CO2 = equilibriumConstantCO2_(temperature, pg, highTemp, spycherPruess2005); // equilibrium constant for CO2 at 1 bar
760 Evaluation phi_CO2 = fugacityCoefficientCO2(temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005); // fugacity coefficient of CO2 for the water-CO2 system
761 Evaluation gammaCO2 = activityCoefficientCO2_(temperature, xCO2, highTemp);
762
763 // In the intermediate temperature range 99-109 C, equilibrium constants and fugacity coeff. are linearly
764 // weighted
765 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
766 const Evaluation weight = (382.15 - temperature) / 10.;
767 const Evaluation& k0_CO2_low = equilibriumConstantCO2_(temperature, pg, false, spycherPruess2005);
768 const Evaluation& phi_CO2_low = fugacityCoefficientCO2(temperature, pg, Evaluation(0.0), false, extrapolate);
769 k0_CO2 = k0_CO2 * (1 - weight) + k0_CO2_low * weight;
770 phi_CO2 = phi_CO2 * (1 - weight) + phi_CO2_low * weight;
771 }
772
773 // Eq. (11)
774 const Evaluation& pg_bar = pg / 1.e5;
775 const Scalar R = IdealGas::R * 10;
776 return phi_CO2 * pg_bar / (55.508 * k0_CO2 * gammaCO2) * exp(-deltaP * v_av_CO2 / (R * temperature));
777 }
778
782 template <class Evaluation>
783 static Evaluation activityCoefficientSalt_(const Evaluation& temperature,
784 const Evaluation& pg,
785 const Evaluation& m_NaCl,
786 const Evaluation& xCO2,
787 const int& activityModel)
788 {
789 OPM_TIMEFUNCTION_LOCAL();
790 // Lambda and xi parameter for either Rumpf et al (1994) (activityModel = 1) or Duan-Sun as modified by Spycher
791 // & Pruess (2009) (activityModel = 2) or Duan & Sun (2003) as given in Spycher & Pruess (2005) (activityModel =
792 // 3)
793 Evaluation lambda;
794 Evaluation xi;
795 Evaluation convTerm;
796 if (activityModel == 1) {
797 lambda = computeLambdaRumpfetal_(temperature);
798 xi = -0.0028 * 3.0;
799 Evaluation m_CO2 = xCO2 * (2 * m_NaCl + 55.508) / (1 - xCO2);
800 convTerm = (1 + (m_CO2 + 2 * m_NaCl) / 55.508) / (1 + m_CO2 / 55.508);
801 }
802 else if (activityModel == 2) {
803 lambda = computeLambdaSpycherPruess2009_(temperature);
804 xi = computeXiSpycherPruess2009_(temperature);
805 convTerm = 1 + 2 * m_NaCl / 55.508;
806 }
807 else if (activityModel == 3) {
808 lambda = computeLambdaDuanSun_(temperature, pg);
809 xi = computeXiDuanSun_(temperature, pg);
810 convTerm = 1.0;
811 }
812 else {
813 throw std::runtime_error("Activity model for salt-out effect has not been implemented!");
814 }
815
816 // Eq. (18)
817 const Evaluation& lnGamma = 2 * lambda * m_NaCl + xi * m_NaCl * m_NaCl;
818
819 // Eq. (18), return activity coeff. on mole-fraction scale
820 return convTerm * exp(lnGamma);
821 }
822
826 template <class Evaluation>
827 static Evaluation computeLambdaSpycherPruess2009_(const Evaluation& temperature)
828 {
829 // Table 1
830 static const Scalar c[3] = { 2.217e-4, 1.074, 2648. };
831
832 // Eq. (19)
833 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
834 }
835
839 template <class Evaluation>
840 static Evaluation computeXiSpycherPruess2009_(const Evaluation& temperature)
841 {
842 // Table 1
843 static const Scalar c[3] = { 1.3e-5, -20.12, 5259. };
844
845 // Eq. (19)
846 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
847 }
848
852 template <class Evaluation>
853 static Evaluation computeLambdaRumpfetal_(const Evaluation& temperature)
854 {
855 // B^(0) below Eq. (A-6)
856 static const Scalar c[4] = { 0.254, -76.82, -10656, 6312e3 };
857
858 return c[0] + c[1] / temperature + c[2] / (temperature * temperature) +
859 c[3] / (temperature * temperature * temperature);
860 }
861
869 template <class Evaluation>
870 static Evaluation computeLambdaDuanSun_(const Evaluation& temperature, const Evaluation& pg)
871 {
872 static const Scalar c[6] =
873 { -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
874
875 Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
876 return c[0] + c[1]*temperature + c[2]/temperature + c[3]*pg_bar/temperature + c[4]*pg_bar/(630.0 - temperature)
877 + c[5]*temperature*log(pg_bar);
878 }
879
887 template <class Evaluation>
888 static Evaluation computeXiDuanSun_(const Evaluation& temperature, const Evaluation& pg)
889 {
890 static const Scalar c[4] =
891 { 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
892
893 Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
894 return c[0] + c[1]*temperature + c[2]*pg_bar/temperature + c[3]*pg_bar/(630.0 - temperature);
895 }
896
903 template <class Evaluation>
904 static Evaluation equilibriumConstantCO2_(const Evaluation& temperature,
905 const Evaluation& pg,
906 const bool& highTemp,
907 bool spycherPruess2005 = false)
908 {
909 OPM_TIMEFUNCTION_LOCAL();
910 Evaluation temperatureCelcius = temperature - 273.15;
911 std::array<Scalar, 4> c;
912 if (highTemp) {
913 c = { 1.668, 3.992e-3, -1.156e-5, 1.593e-9 };
914 }
915 else {
916 // For temperature below 31 C and pressures above saturation pressure, separate parameters are needed
917 Evaluation psat = CO2::vaporPressure(temperature);
918 if (temperatureCelcius < 31 && pg > psat && !spycherPruess2005) {
919 c = { 1.169, 1.368e-2, -5.38e-5, 0.0 };
920 }
921 else {
922 c = { 1.189, 1.304e-2, -5.446e-5, 0.0 };
923 }
924 }
925 Evaluation logk0_CO2 = c[0] + temperatureCelcius * (c[1] + temperatureCelcius *
926 (c[2] + temperatureCelcius * c[3]));
927 Evaluation k0_CO2 = pow(10.0, logk0_CO2);
928 return k0_CO2;
929 }
930
937 template <class Evaluation>
938 static Evaluation equilibriumConstantH2O_(const Evaluation& temperature, const bool& highTemp)
939 {
940 Evaluation temperatureCelcius = temperature - 273.15;
941 std::array<Scalar, 5> c;
942 if (highTemp){
943 c = { -2.1077, 2.8127e-2, -8.4298e-5, 1.4969e-7, -1.1812e-10 };
944 }
945 else {
946 c = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7, 0.0 };
947 }
948 Evaluation logk0_H2O = c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
949 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
950 return pow(10.0, logk0_H2O);
951 }
952
953};
954
955} // namespace BinaryCoeff
956} // namespace Opm
957
958#endif
Relations valid for an ideal gas.
Some templates to wrap the valgrind client request macros.
Binary coefficients for brine and CO2.
Definition Brine_CO2.hpp:45
static void calculateMoleFractions(const Evaluation &temperature, const Evaluation &pg, const Evaluation &salinity, const int knownPhaseIdx, Evaluation &xlCO2, Evaluation &ygH2O, const int &activityModel, bool extrapolate=false)
Returns the mol (!) fraction of CO2 in the liquid phase and the mol_ (!) fraction of H2O in the gas p...
Definition Brine_CO2.hpp:100
static Evaluation fugacityCoefficientH2O(const Evaluation &temperature, const Evaluation &pg, const Evaluation &yH2O, const bool highTemp, bool extrapolate=false, bool spycherPruess2005=false)
Returns the fugacity coefficient of the H2O component in a water-CO2 mixture.
Definition Brine_CO2.hpp:240
static Evaluation gasDiffCoeff(const Evaluation &temperature, const Evaluation &pressure, bool extrapolate=false)
Binary diffusion coefficent [m^2/s] of water in the CO2 phase.
Definition Brine_CO2.hpp:59
static Evaluation fugacityCoefficientCO2(const Evaluation &temperature, const Evaluation &pg, const Evaluation &yH2O, const bool highTemp, bool extrapolate=false, bool spycherPruess2005=false)
Returns the fugacity coefficient of the CO2 component in a water-CO2 mixture.
Definition Brine_CO2.hpp:185
static Evaluation liquidDiffCoeff(const Evaluation &, const Evaluation &)
Binary diffusion coefficent [m^2/s] of CO2 in the brine phase.
Definition Brine_CO2.hpp:76
static Evaluation henry(const Evaluation &temperature, bool extrapolate=false)
Henry coefficent for CO2 in brine.
Definition Brine_CO2.hpp:173
static Evaluation gasViscosity(Evaluation temperature, const Evaluation &pressure, bool extrapolate=false)
The dynamic viscosity [Pa s] of CO2.
Definition CO2.hpp:208
static Evaluation vaporPressure(const Evaluation &T)
Returns the pressure [Pa] at CO2's triple point.
Definition CO2.hpp:135
static Scalar molarMass()
The mass in [kg] of one mole of CO2.
Definition CO2.hpp:71
static Evaluation gasDensity(const Evaluation &temperature, const Evaluation &pressure, bool extrapolate=false)
The density of CO2 at a given pressure and temperature [kg/m^3].
Definition CO2.hpp:194
static Evaluation vaporPressure(Evaluation temperature)
The vapor pressure in of pure water at a given temperature.
Definition H2O.hpp:143
static const Scalar molarMass()
The molar mass in of water.
Definition H2O.hpp:85
Relations valid for an ideal gas.
Definition IdealGas.hpp:38
static const Scalar R
The ideal gas constant .
Definition IdealGas.hpp:41
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30