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#include <opm/common/utility/gpuDecorators.hpp>
35
36#include <array>
37
38namespace Opm {
39namespace BinaryCoeff {
40
45template<class Scalar, class H2O, class CO2, bool verbose = true>
46class Brine_CO2 {
47 typedef ::Opm::IdealGas<Scalar> IdealGas;
48 static const int liquidPhaseIdx = 0; // index of the liquid phase
49 static const int gasPhaseIdx = 1; // index of the gas phase
50
51public:
59 template <class Evaluation, class CO2Params>
60 OPM_HOST_DEVICE static Evaluation gasDiffCoeff(const CO2Params& params, const Evaluation& temperature, const Evaluation& pressure, bool extrapolate = false)
61 {
62 //Diffusion coefficient of water in the CO2 phase
63 Scalar k = 1.3806504e-23; // Boltzmann constant
64 Scalar c = 4; // slip parameter, can vary between 4 (slip condition) and 6 (stick condition)
65 Scalar R_h = 1.72e-10; // hydrodynamic radius of the solute
66 const Evaluation& mu = CO2::gasViscosity(params, temperature, pressure, extrapolate); // CO2 viscosity
67 return k / (c * M_PI * R_h) * (temperature / mu);
68 }
69
76 template <class Evaluation>
77 OPM_HOST_DEVICE static Evaluation liquidDiffCoeff(const Evaluation& /*temperature*/, const Evaluation& /*pressure*/)
78 {
79 //Diffusion coefficient of CO2 in the brine phase
80 return 2e-9;
81 }
82
100 template <class Evaluation, class CO2Params>
101 OPM_HOST_DEVICE static void calculateMoleFractions(const CO2Params& params,
102 const Evaluation& temperature,
103 const Evaluation& pg,
104 const Evaluation& salinity,
105 const int knownPhaseIdx,
106 Evaluation& xlCO2,
107 Evaluation& ygH2O,
108 const int& activityModel,
109 bool extrapolate = false)
110 {
111 OPM_TIMEFUNCTION_LOCAL();
112
113 // Iterate or not?
114 bool iterate = false;
115 if ((activityModel == 1 && salinity > 0.0) || (activityModel == 2 && temperature > 372.15)) {
116 iterate = true;
117 }
118
119 // If both phases are present the mole fractions in each phase can be calculate with the mutual solubility
120 // function
121 if (knownPhaseIdx < 0) {
122 Evaluation molalityNaCl = massFracToMolality_(salinity); // mass fraction to molality of NaCl
123
124 // Duan-Sun model as given in Spycher & Pruess (2005) have a different fugacity coefficient formula and
125 // activity coefficient definition (not a true activity coefficient but a ratio).
126 // Technically only valid below T = 100 C, but we use low-temp. parameters and formulas even above 100 C as
127 // an approximation.
128 if (activityModel == 3) {
129 auto [xCO2, yH2O] = mutualSolubilitySpycherPruess2005_(params, temperature, pg, molalityNaCl, extrapolate);
130 xlCO2 = xCO2;
131 ygH2O = yH2O;
132
133 }
134 else {
135 // Fixed-point iterations to calculate solubility
136 if (iterate) {
137 auto [xCO2, yH2O] = fixPointIterSolubility_(params, temperature, pg, molalityNaCl, activityModel, extrapolate);
138 xlCO2 = xCO2;
139 ygH2O = yH2O;
140 }
141
142 // Solve mutual solubility equation with back substitution (no need for iterations)
143 else {
144 auto [xCO2, yH2O] = nonIterSolubility_(params, temperature, pg, molalityNaCl, activityModel, extrapolate);
145 xlCO2 = xCO2;
146 ygH2O = yH2O;
147 }
148 }
149 }
150
151 // if only liquid phase is present the mole fraction of CO2 in brine is given and
152 // and the virtual equilibrium mole fraction of water in the non-existing gas phase can be estimated
153 // with the mutual solubility function
154 else if (knownPhaseIdx == liquidPhaseIdx && activityModel == 3) {
155 Evaluation x_NaCl = salinityToMolFrac_(salinity);
156 const Evaluation& A = computeA_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
157 ygH2O = A * (1 - xlCO2 - x_NaCl);
158 }
159
160 // if only gas phase is present the mole fraction of water in the gas phase is given and
161 // and the virtual equilibrium mole fraction of CO2 in the non-existing liquid phase can be estimated
162 // with the mutual solubility function
163 else if (knownPhaseIdx == gasPhaseIdx && activityModel == 3) {
164 //y_H2o = fluidstate.
165 Evaluation x_NaCl = salinityToMolFrac_(salinity);
166 const Evaluation& A = computeA_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
167 xlCO2 = 1 - x_NaCl - ygH2O / A;
168 }
169 }
170
174 template <class Evaluation>
175 static Evaluation henry(const Evaluation& temperature, bool extrapolate = false)
176 { return fugacityCoefficientCO2(temperature, /*pressure=*/1e5, extrapolate)*1e5; }
177
186 template <class Evaluation, class CO2Params>
187 static Evaluation fugacityCoefficientCO2(const CO2Params& params,
188 const Evaluation& temperature,
189 const Evaluation& pg,
190 const Evaluation& yH2O,
191 const bool highTemp,
192 bool extrapolate = false,
193 bool spycherPruess2005 = false)
194 {
195 OPM_TIMEFUNCTION_LOCAL();
196 Valgrind::CheckDefined(temperature);
197 Valgrind::CheckDefined(pg);
198
199 const Evaluation V = 1 / (CO2::gasDensity(params, temperature, pg, extrapolate) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
200 const Evaluation pg_bar = pg / 1.e5; // gas phase pressure in bar
201 const Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
202
203 // Parameters in Redlich-Kwong equation
204 const Evaluation a_CO2 = aCO2_(temperature, highTemp);
205 const Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
206 const Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
207 const Scalar b_CO2 = bCO2_(highTemp);
208 const Evaluation b_mix = bMix_(yH2O, highTemp);
209 const Evaluation Rt15 = R * pow(temperature, 1.5);
210
211 Evaluation lnPhiCO2;
212 if (spycherPruess2005) {
213 const Evaluation logVpb_V = log((V + b_CO2) / V);
214 lnPhiCO2 = log(V / (V - b_CO2));
215 lnPhiCO2 += b_CO2 / (V - b_CO2);
216 lnPhiCO2 -= 2 * a_CO2 / (Rt15 * b_CO2) * logVpb_V;
217 lnPhiCO2 +=
218 a_CO2 * b_CO2
219 / (Rt15
220 * b_CO2
221 * b_CO2)
222 * (logVpb_V
223 - b_CO2 / (V + b_CO2));
224 lnPhiCO2 -= log(pg_bar * V / (R * temperature));
225 }
226 else {
227 lnPhiCO2 = (b_CO2 / b_mix) * (pg_bar * V / (R * temperature) - 1);
228 lnPhiCO2 -= log(pg_bar * (V - b_mix) / (R * temperature));
229 lnPhiCO2 += (2 * (yH2O * a_CO2_H2O + (1 - yH2O) * a_CO2) / a_mix - (b_CO2 / b_mix)) *
230 a_mix / (b_mix * Rt15) * log(V / (V + b_mix));
231 }
232 return exp(lnPhiCO2); // fugacity coefficient of CO2
233 }
234
243 template <class Evaluation, class CO2Params>
244 static Evaluation fugacityCoefficientH2O(const CO2Params& params,
245 const Evaluation& temperature,
246 const Evaluation& pg,
247 const Evaluation& yH2O,
248 const bool highTemp,
249 bool extrapolate = false,
250 bool spycherPruess2005 = false)
251 {
252 OPM_TIMEFUNCTION_LOCAL();
253 Valgrind::CheckDefined(temperature);
254 Valgrind::CheckDefined(pg);
255
256 const Evaluation& V = 1 / (CO2::gasDensity(params, temperature, pg, extrapolate) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
257 const Evaluation& pg_bar = pg / 1.e5; // gas phase pressure in bar
258 const Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
259
260 // Mixture parameter of Redlich-Kwong equation
261 const Evaluation a_H2O = aH2O_(temperature, highTemp);
262 const Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
263 const Evaluation a_mix = aMix_(temperature, yH2O, highTemp);
264 const Scalar b_H2O = bH2O_(highTemp);
265 const Evaluation b_mix = bMix_(yH2O, highTemp);
266 const Evaluation Rt15 = R * pow(temperature, 1.5);
267
268 Evaluation lnPhiH2O;
269 if (spycherPruess2005) {
270 const Evaluation logVpb_V = log((V + b_mix) / V);
271 lnPhiH2O =
272 log(V/(V - b_mix))
273 + b_H2O/(V - b_mix) - 2*a_CO2_H2O
274 / (Rt15*b_mix)*logVpb_V
275 + a_mix*b_H2O/(Rt15*b_mix*b_mix)
276 *(logVpb_V - b_mix/(V + b_mix))
277 - log(pg_bar*V/(R*temperature));
278 }
279 else {
280 lnPhiH2O = (b_H2O / b_mix) * (pg_bar * V / (R * temperature) - 1);
281 lnPhiH2O -= log(pg_bar * (V - b_mix) / (R * temperature));
282 lnPhiH2O += (2 * (yH2O * a_H2O + (1 - yH2O) * a_CO2_H2O) / a_mix - (b_H2O / b_mix)) *
283 a_mix / (b_mix * Rt15) * log(V / (V + b_mix));
284 }
285 return exp(lnPhiH2O); // fugacity coefficient of H2O
286 }
287
288private:
292 template <class Evaluation>
293 OPM_HOST_DEVICE static Evaluation aCO2_(const Evaluation& temperature, const bool& highTemp)
294 {
295 if (highTemp) {
296 return 8.008e7 - 4.984e4 * temperature;
297 }
298 else {
299 return 7.54e7 - 4.13e4 * temperature;
300 }
301 }
302
306 template <class Evaluation>
307 OPM_HOST_DEVICE static Evaluation aH2O_(const Evaluation& temperature, const bool& highTemp)
308 {
309 if (highTemp) {
310 return 1.337e8 - 1.4e4 * temperature;
311 }
312 else {
313 return 0.0;
314 }
315 }
316
320 template <class Evaluation>
321 OPM_HOST_DEVICE static Evaluation aCO2_H2O_(const Evaluation& temperature, const Evaluation& yH2O, const bool& highTemp)
322 {
323 if (highTemp) {
324 // Pure parameters
325 Evaluation aCO2 = aCO2_(temperature, highTemp);
326 Evaluation aH2O = aH2O_(temperature, highTemp);
327
328 // Mixture Eq. (A-6)
329 Evaluation K_CO2_H2O = 0.4228 - 7.422e-4 * temperature;
330 Evaluation K_H2O_CO2 = 1.427e-2 - 4.037e-4 * temperature;
331 Evaluation k_CO2_H2O = yH2O * K_H2O_CO2 + (1 - yH2O) * K_CO2_H2O;
332
333 // Eq. (A-5)
334 return sqrt(aCO2 * aH2O) * (1 - k_CO2_H2O);
335 }
336 else {
337 return 7.89e7;
338 }
339 }
340
344 template <class Evaluation>
345 OPM_HOST_DEVICE static Evaluation aMix_(const Evaluation& temperature, const Evaluation& yH2O, const bool& highTemp)
346 {
347 if (highTemp) {
348 // Parameters
349 Evaluation aCO2 = aCO2_(temperature, highTemp);
350 Evaluation aH2O = aH2O_(temperature, highTemp);
351 Evaluation a_CO2_H2O = aCO2_H2O_(temperature, yH2O, highTemp);
352
353 return yH2O * yH2O * aH2O + 2 * yH2O * (1 - yH2O) * a_CO2_H2O + (1 - yH2O) * (1 - yH2O) * aCO2;
354 }
355 else {
356 return aCO2_(temperature, highTemp);
357 }
358 }
359
363 OPM_HOST_DEVICE static Scalar bCO2_(const bool& highTemp)
364 {
365 if (highTemp) {
366 return 28.25;
367 }
368 else {
369 return 27.8;
370 }
371 }
372
376 OPM_HOST_DEVICE static Scalar bH2O_(const bool& highTemp)
377 {
378 if (highTemp) {
379 return 15.7;
380 }
381 else {
382 return 18.18;
383 }
384 }
385
389 template <class Evaluation>
390 OPM_HOST_DEVICE static Evaluation bMix_(const Evaluation& yH2O, const bool& highTemp)
391 {
392 if (highTemp) {
393 // Parameters
394 Scalar bCO2 = bCO2_(highTemp);
395 Scalar bH2O = bH2O_(highTemp);
396
397 return yH2O * bH2O + (1 - yH2O) * bCO2;
398 }
399 else {
400 return bCO2_(highTemp);
401 }
402 }
403
407 template <class Evaluation>
408 OPM_HOST_DEVICE static Evaluation V_avg_CO2_(const Evaluation& temperature, const bool& highTemp)
409 {
410 if (highTemp && (temperature > 373.15)) {
411 return 32.6 + 3.413e-2 * (temperature - 373.15);
412 }
413 else {
414 return 32.6;
415 }
416 }
417
421 template <class Evaluation>
422 OPM_HOST_DEVICE static Evaluation V_avg_H2O_(const Evaluation& temperature, const bool& highTemp)
423 {
424 if (highTemp && (temperature > 373.15)) {
425 return 18.1 + 3.137e-2 * (temperature - 373.15);
426 }
427 else {
428 return 18.1;
429 }
430 }
431
435 template <class Evaluation>
436 OPM_HOST_DEVICE static Evaluation AM_(const Evaluation& temperature, const bool& highTemp)
437 {
438 if (highTemp && temperature > 373.15) {
439 Evaluation deltaTk = temperature - 373.15;
440 return deltaTk * (-3.084e-2 + 1.927e-5 * deltaTk);
441 }
442 else {
443 return 0.0;
444 }
445 }
446
450 template <class Evaluation>
451 OPM_HOST_DEVICE static Evaluation Pref_(const Evaluation& temperature, const bool& highTemp)
452 {
453 if (highTemp && temperature > 373.15) {
454 const Evaluation& temperatureCelcius = temperature - 273.15;
455 static const Scalar c[5] = { -1.9906e-1, 2.0471e-3, 1.0152e-4, -1.4234e-6, 1.4168e-8 };
456 return c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
457 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
458 }
459 else {
460 return 1.0;
461 }
462 }
463
467 template <class Evaluation>
468 OPM_HOST_DEVICE static Evaluation activityCoefficientCO2_(const Evaluation& temperature,
469 const Evaluation& xCO2,
470 const bool& highTemp)
471 {
472 if (highTemp) {
473 // Eq. (13)
474 Evaluation AM = AM_(temperature, highTemp);
475 Evaluation lnGammaCO2 = 2 * AM * xCO2 * (1 - xCO2) * (1 - xCO2);
476 return exp(lnGammaCO2);
477 }
478 else {
479 return 1.0;
480 }
481 }
482
486 template <class Evaluation>
487 OPM_HOST_DEVICE static Evaluation activityCoefficientH2O_(const Evaluation& temperature,
488 const Evaluation& xCO2,
489 const bool& highTemp)
490 {
491 if (highTemp) {
492 // Eq. (12)
493 Evaluation AM = AM_(temperature, highTemp);
494 Evaluation lnGammaH2O = (1 - 2 * (1 - xCO2)) * AM * xCO2 * xCO2;
495 return exp(lnGammaH2O);
496 }
497 else {
498 return 1.0;
499 }
500 }
501
507 template <class Evaluation>
508 OPM_HOST_DEVICE static Evaluation salinityToMolFrac_(const Evaluation& salinity) {
509 OPM_TIMEFUNCTION_LOCAL();
510 const Scalar Mw = H2O::molarMass(); /* molecular weight of water [kg/mol] */
511 const Scalar Ms = 58.44e-3; /* molecular weight of NaCl [kg/mol] */
512
513 const Evaluation X_NaCl = salinity;
514 /* salinity: conversion from mass fraction to mol fraction */
515 const Evaluation x_NaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
516 return x_NaCl;
517 }
518
524#if 0
525 template <class Evaluation>
526 OPM_HOST_DEVICE static Evaluation moleFracToMolality_(const Evaluation& x_NaCl)
527 {
528 // conversion from mol fraction to molality (dissolved CO2 neglected)
529 return 55.508 * x_NaCl / (1 - x_NaCl);
530 }
531#endif
532
533 template <class Evaluation>
534 OPM_HOST_DEVICE static Evaluation massFracToMolality_(const Evaluation& X_NaCl)
535 {
536 const Scalar MmNaCl = 58.44e-3;
537 return X_NaCl / (MmNaCl * (1 - X_NaCl));
538 }
539
545 template <class Evaluation>
546 OPM_HOST_DEVICE static Evaluation molalityToMoleFrac_(const Evaluation& m_NaCl)
547 {
548 // conversion from molality to mole fractio (dissolved CO2 neglected)
549 return m_NaCl / (55.508 + m_NaCl);
550 }
551
555 template <class Evaluation, class CO2Parameters>
556 OPM_HOST_DEVICE static std::pair<Evaluation, Evaluation> fixPointIterSolubility_(const CO2Parameters& params,
557 const Evaluation& temperature,
558 const Evaluation& pg,
559 const Evaluation& m_NaCl,
560 const int& activityModel,
561 bool extrapolate = false)
562 {
563 OPM_TIMEFUNCTION_LOCAL();
564 // Start point for fixed-point iterations as recommended below in section 2.2
565 Evaluation yH2O = H2O::vaporPressure(temperature) / pg; // ideal mixing
566 Evaluation xCO2 = 0.009; // same as ~0.5 mol/kg
567 Evaluation gammaNaCl = 1.0; // default salt activity coeff = 1.0
568
569 // We can pre-calculate Duan-Sun, Spycher & Pruess (2009) salt activity coeff.
570 if (m_NaCl > 0.0 && activityModel == 2) {
571 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
572 }
573
574 // Options
575 int max_iter = 100;
576 Scalar tol = 1e-8;
577 bool highTemp = true;
578 if (activityModel == 1) {
579 highTemp = false;
580 }
581 const bool iterate = true;
582
583 // Fixed-point loop x_i+1 = F(x_i)
584 for (int i = 0; i < max_iter; ++i) {
585 // Calculate activity coefficient for Rumpf et al (1994) model
586 if (m_NaCl > 0.0 && activityModel == 1) {
587 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, xCO2, activityModel);
588 }
589
590 // F(x_i) is the mutual solubilities
591 auto [xCO2_new, yH2O_new] = mutualSolubility_(params, temperature, pg, xCO2, yH2O, m_NaCl, gammaNaCl, highTemp,
592 iterate, extrapolate);
593
594 // Check for convergence
595 if (abs(xCO2_new - xCO2) < tol && abs(yH2O_new - yH2O) < tol) {
596 xCO2 = xCO2_new;
597 yH2O = yH2O_new;
598 break;
599 }
600
601 // Else update mole fractions for next iteration
602 else {
603 xCO2 = xCO2_new;
604 yH2O = yH2O_new;
605 }
606 }
607
608 return {xCO2, yH2O};
609 }
610
614 template <class Evaluation, class CO2Parameters>
615 OPM_HOST_DEVICE static std::pair<Evaluation, Evaluation> nonIterSolubility_(const CO2Parameters& params,
616 const Evaluation& temperature,
617 const Evaluation& pg,
618 const Evaluation& m_NaCl,
619 const int& activityModel,
620 bool extrapolate = false)
621 {
622 // Calculate activity coefficient for salt
623 Evaluation gammaNaCl = 1.0;
624 if (m_NaCl > 0.0 && activityModel > 0 && activityModel < 3) {
625 gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), activityModel);
626 }
627
628 // Calculate mutual solubility.
629 // Note that we don't use xCO2 and yH2O input in low-temperature case, so we set them to 0.0
630 const bool highTemp = false;
631 const bool iterate = false;
632 auto [xCO2, yH2O] = mutualSolubility_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0), m_NaCl, gammaNaCl,
633 highTemp, iterate, extrapolate);
634
635 return {xCO2, yH2O};
636 }
637
641 template <class Evaluation, class CO2Parameters>
642 OPM_HOST_DEVICE static std::pair<Evaluation, Evaluation> mutualSolubility_(const CO2Parameters& params,
643 const Evaluation& temperature,
644 const Evaluation& pg,
645 const Evaluation& xCO2,
646 const Evaluation& yH2O,
647 const Evaluation& m_NaCl,
648 const Evaluation& gammaNaCl,
649 const bool& highTemp,
650 const bool& iterate,
651 bool extrapolate = false)
652 {
653 // Calculate A and B (without salt effect); Eqs. (8) and (9)
654 const Evaluation& A = computeA_(params, temperature, pg, yH2O, xCO2, highTemp, extrapolate);
655 Evaluation B = computeB_(params, temperature, pg, yH2O, xCO2, highTemp, extrapolate);
656
657 // Add salt effect to B, Eq. (17)
658 B /= gammaNaCl;
659
660 // Compute yH2O and xCO2, Eqs. (B-7) and (B-2)
661 Evaluation yH2O_new = (1. - B) * 55.508 / ((1. / A - B) * (2 * m_NaCl + 55.508) + 2 * m_NaCl * B);
662 Evaluation xCO2_new;
663 if (iterate) {
664 xCO2_new = B * (1 - yH2O);
665 }
666 else {
667 xCO2_new = B * (1 - yH2O_new);
668 }
669
670 return {xCO2_new, yH2O_new};
671 }
672
676 template <class Evaluation, class CO2Parameters>
677 OPM_HOST_DEVICE static std::pair<Evaluation, Evaluation> mutualSolubilitySpycherPruess2005_(const CO2Parameters& params,
678 const Evaluation& temperature,
679 const Evaluation& pg,
680 const Evaluation& m_NaCl,
681 bool extrapolate = false)
682 {
683 // Calculate A and B (without salt effect); Eqs. (8) and (9)
684 const Evaluation& A = computeA_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
685 const Evaluation& B = computeB_(params, temperature, pg, Evaluation(0.0), Evaluation(0.0), false, extrapolate, true);
686
687 // Mole fractions and molality in pure water
688 Evaluation yH2O = (1 - B) / (1. / A - B);
689 Evaluation xCO2 = B * (1 - yH2O);
690
691 // Modifiy mole fractions with Duan-Sun "activity coefficient" if salt is involved
692 if (m_NaCl > 0.0) {
693 const Evaluation& gammaNaCl = activityCoefficientSalt_(temperature, pg, m_NaCl, Evaluation(0.0), 3);
694
695 // Molality with salt
696 Evaluation mCO2 = (xCO2 * 55.508) / (1 - xCO2); // pure water
697 mCO2 /= gammaNaCl;
698 xCO2 = mCO2 / (m_NaCl + 55.508 + mCO2);
699
700 // new yH2O with salt
701 const Evaluation& xNaCl = molalityToMoleFrac_(m_NaCl);
702 yH2O = A * (1 - xCO2 - xNaCl);
703 }
704
705 return {xCO2, yH2O};
706 }
707
716 template <class Evaluation, class CO2Params>
717 OPM_HOST_DEVICE static Evaluation computeA_(const CO2Params& params,
718 const Evaluation& temperature,
719 const Evaluation& pg,
720 const Evaluation& yH2O,
721 const Evaluation& xCO2,
722 const bool& highTemp,
723 bool extrapolate = false,
724 bool spycherPruess2005 = false)
725 {
726 OPM_TIMEFUNCTION_LOCAL();
727 // Intermediate calculations
728 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp); // pressure range [bar] from pref to pg[bar]
729 Evaluation v_av_H2O = V_avg_H2O_(temperature, highTemp); // average partial molar volume of H2O [cm^3/mol]
730 Evaluation k0_H2O = equilibriumConstantH2O_(temperature, highTemp); // equilibrium constant for H2O at 1 bar
731 Evaluation phi_H2O = fugacityCoefficientH2O(params, temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005); // fugacity coefficient of H2O for the water-CO2 system
732 Evaluation gammaH2O = activityCoefficientH2O_(temperature, xCO2, highTemp);
733
734 // In the intermediate temperature range 99-109 C, equilibrium constants and fugacity coeff. are linearly
735 // weighted
736 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
737 const Evaluation weight = (382.15 - temperature) / 10.;
738 const Evaluation& k0_H2O_low = equilibriumConstantH2O_(temperature, false);
739 const Evaluation& phi_H2O_low = fugacityCoefficientH2O(params, temperature, pg, Evaluation(0.0), false, extrapolate);
740 k0_H2O = k0_H2O * (1 - weight) + k0_H2O_low * weight;
741 phi_H2O = phi_H2O * (1 - weight) + phi_H2O_low * weight;
742 }
743
744 // Eq. (10)
745 const Evaluation& pg_bar = pg / 1.e5;
746 Scalar R = IdealGas::R * 10;
747 return k0_H2O * gammaH2O / (phi_H2O * pg_bar) * exp(deltaP * v_av_H2O / (R * temperature));
748 }
749
758 template <class Evaluation, class CO2Parameters>
759 static Evaluation computeB_(const CO2Parameters& params,
760 const Evaluation& temperature,
761 const Evaluation& pg,
762 const Evaluation& yH2O,
763 const Evaluation& xCO2,
764 const bool& highTemp,
765 bool extrapolate = false,
766 bool spycherPruess2005 = false)
767 {
768 OPM_TIMEFUNCTION_LOCAL();
769 // Intermediate calculations
770 const Evaluation& deltaP = pg / 1e5 - Pref_(temperature, highTemp); // pressure range [bar] from pref to pg[bar]
771 Evaluation v_av_CO2 = V_avg_CO2_(temperature, highTemp); // average partial molar volume of CO2 [cm^3/mol]
772 Evaluation k0_CO2 = equilibriumConstantCO2_(temperature, pg, highTemp, spycherPruess2005); // equilibrium constant for CO2 at 1 bar
773 Evaluation phi_CO2 = fugacityCoefficientCO2(params, temperature, pg, yH2O, highTemp, extrapolate, spycherPruess2005); // fugacity coefficient of CO2 for the water-CO2 system
774 Evaluation gammaCO2 = activityCoefficientCO2_(temperature, xCO2, highTemp);
775
776 // In the intermediate temperature range 99-109 C, equilibrium constants and fugacity coeff. are linearly
777 // weighted
778 if ( temperature > 372.15 && temperature < 382.15 && !spycherPruess2005) {
779 const Evaluation weight = (382.15 - temperature) / 10.;
780 const Evaluation& k0_CO2_low = equilibriumConstantCO2_(temperature, pg, false, spycherPruess2005);
781 const Evaluation& phi_CO2_low = fugacityCoefficientCO2(params, temperature, pg, Evaluation(0.0), false, extrapolate);
782 k0_CO2 = k0_CO2 * (1 - weight) + k0_CO2_low * weight;
783 phi_CO2 = phi_CO2 * (1 - weight) + phi_CO2_low * weight;
784 }
785
786 // Eq. (11)
787 const Evaluation& pg_bar = pg / 1.e5;
788 const Scalar R = IdealGas::R * 10;
789 return phi_CO2 * pg_bar / (55.508 * k0_CO2 * gammaCO2) * exp(-deltaP * v_av_CO2 / (R * temperature));
790 }
791
795 template <class Evaluation>
796 OPM_HOST_DEVICE static Evaluation activityCoefficientSalt_(const Evaluation& temperature,
797 const Evaluation& pg,
798 const Evaluation& m_NaCl,
799 const Evaluation& xCO2,
800 const int& activityModel)
801 {
802 OPM_TIMEFUNCTION_LOCAL();
803 // Lambda and xi parameter for either Rumpf et al (1994) (activityModel = 1) or Duan-Sun as modified by Spycher
804 // & Pruess (2009) (activityModel = 2) or Duan & Sun (2003) as given in Spycher & Pruess (2005) (activityModel =
805 // 3)
806 Evaluation lambda;
807 Evaluation xi;
808 Evaluation convTerm;
809 if (activityModel == 1) {
810 lambda = computeLambdaRumpfetal_(temperature);
811 xi = -0.0028 * 3.0;
812 Evaluation m_CO2 = xCO2 * (2 * m_NaCl + 55.508) / (1 - xCO2);
813 convTerm = (1 + (m_CO2 + 2 * m_NaCl) / 55.508) / (1 + m_CO2 / 55.508);
814 }
815 else if (activityModel == 2) {
816 lambda = computeLambdaSpycherPruess2009_(temperature);
817 xi = computeXiSpycherPruess2009_(temperature);
818 convTerm = 1 + 2 * m_NaCl / 55.508;
819 }
820 else if (activityModel == 3) {
821 lambda = computeLambdaDuanSun_(temperature, pg);
822 xi = computeXiDuanSun_(temperature, pg);
823 convTerm = 1.0;
824 }
825 else {
826 throw std::runtime_error("Activity model for salt-out effect has not been implemented!");
827 }
828
829 // Eq. (18)
830 const Evaluation& lnGamma = 2 * lambda * m_NaCl + xi * m_NaCl * m_NaCl;
831
832 // Eq. (18), return activity coeff. on mole-fraction scale
833 return convTerm * exp(lnGamma);
834 }
835
839 template <class Evaluation>
840 OPM_HOST_DEVICE static Evaluation computeLambdaSpycherPruess2009_(const Evaluation& temperature)
841 {
842 // Table 1
843 static const Scalar c[3] = { 2.217e-4, 1.074, 2648. };
844
845 // Eq. (19)
846 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
847 }
848
852 template <class Evaluation>
853 OPM_HOST_DEVICE static Evaluation computeXiSpycherPruess2009_(const Evaluation& temperature)
854 {
855 // Table 1
856 static const Scalar c[3] = { 1.3e-5, -20.12, 5259. };
857
858 // Eq. (19)
859 return c[0] * temperature + c[1] / temperature + c[2] / (temperature * temperature);
860 }
861
865 template <class Evaluation>
866 OPM_HOST_DEVICE static Evaluation computeLambdaRumpfetal_(const Evaluation& temperature)
867 {
868 // B^(0) below Eq. (A-6)
869 static const Scalar c[4] = { 0.254, -76.82, -10656, 6312e3 };
870
871 return c[0] + c[1] / temperature + c[2] / (temperature * temperature) +
872 c[3] / (temperature * temperature * temperature);
873 }
874
882 template <class Evaluation>
883 OPM_HOST_DEVICE static Evaluation computeLambdaDuanSun_(const Evaluation& temperature, const Evaluation& pg)
884 {
885 static const Scalar c[6] =
886 { -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
887
888 Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
889 return c[0] + c[1]*temperature + c[2]/temperature + c[3]*pg_bar/temperature + c[4]*pg_bar/(630.0 - temperature)
890 + c[5]*temperature*log(pg_bar);
891 }
892
900 template <class Evaluation>
901 OPM_HOST_DEVICE static Evaluation computeXiDuanSun_(const Evaluation& temperature, const Evaluation& pg)
902 {
903 static const Scalar c[4] =
904 { 3.36389723E-4, -1.98298980E-5, 2.12220830E-3, -5.24873303E-3 };
905
906 Evaluation pg_bar = pg / 1.0E5; /* conversion from Pa to bar */
907 return c[0] + c[1]*temperature + c[2]*pg_bar/temperature + c[3]*pg_bar/(630.0 - temperature);
908 }
909
916 template <class Evaluation>
917 static Evaluation equilibriumConstantCO2_(const Evaluation& temperature,
918 const Evaluation& pg,
919 const bool& highTemp,
920 bool spycherPruess2005 = false)
921 {
922 OPM_TIMEFUNCTION_LOCAL();
923 Evaluation temperatureCelcius = temperature - 273.15;
924 std::array<Scalar, 4> c;
925 if (highTemp) {
926 c = { 1.668, 3.992e-3, -1.156e-5, 1.593e-9 };
927 }
928 else {
929 // For temperature below critical temperature and pressures above saturation pressure, separate parameters are needed
930 bool model1 = temperature < CO2::criticalTemperature() && !spycherPruess2005;
931 if (model1) {
932 // Computing the vapor pressure is not trivial and is also not defined for T > criticalTemperature
933 Evaluation psat = CO2::vaporPressure(temperature);
934 model1 = pg > psat;
935 }
936 if (model1) {
937 c = { 1.169, 1.368e-2, -5.38e-5, 0.0 };
938 }
939 else {
940 c = { 1.189, 1.304e-2, -5.446e-5, 0.0 };
941 }
942 }
943 Evaluation logk0_CO2 = c[0] + temperatureCelcius * (c[1] + temperatureCelcius *
944 (c[2] + temperatureCelcius * c[3]));
945 Evaluation k0_CO2 = pow(10.0, logk0_CO2);
946 return k0_CO2;
947 }
948
955 template <class Evaluation>
956 OPM_HOST_DEVICE static Evaluation equilibriumConstantH2O_(const Evaluation& temperature, const bool& highTemp)
957 {
958 Evaluation temperatureCelcius = temperature - 273.15;
959 std::array<Scalar, 5> c;
960 if (highTemp){
961 c = { -2.1077, 2.8127e-2, -8.4298e-5, 1.4969e-7, -1.1812e-10 };
962 }
963 else {
964 c = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7, 0.0 };
965 }
966 Evaluation logk0_H2O = c[0] + temperatureCelcius * (c[1] + temperatureCelcius * (c[2] +
967 temperatureCelcius * (c[3] + temperatureCelcius * c[4])));
968 return pow(10.0, logk0_H2O);
969 }
970
971};
972
973} // namespace BinaryCoeff
974} // namespace Opm
975
976#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:46
static OPM_HOST_DEVICE Evaluation liquidDiffCoeff(const Evaluation &, const Evaluation &)
Binary diffusion coefficent [m^2/s] of CO2 in the brine phase.
Definition Brine_CO2.hpp:77
static OPM_HOST_DEVICE void calculateMoleFractions(const CO2Params &params, 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:101
static Evaluation fugacityCoefficientH2O(const CO2Params &params, 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:244
static OPM_HOST_DEVICE Evaluation gasDiffCoeff(const CO2Params &params, 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:60
static Evaluation fugacityCoefficientCO2(const CO2Params &params, 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:187
static Evaluation henry(const Evaluation &temperature, bool extrapolate=false)
Henry coefficent for CO2 in brine.
Definition Brine_CO2.hpp:175
static OPM_HOST_DEVICE Evaluation gasDensity(const Params &params, 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:222
static OPM_HOST_DEVICE Scalar molarMass()
The mass in [kg] of one mole of CO2.
Definition CO2.hpp:73
static OPM_HOST_DEVICE Scalar criticalTemperature()
Returns the critical temperature [K] of CO2.
Definition CO2.hpp:79
static OPM_HOST_DEVICE Evaluation vaporPressure(const Evaluation &T)
Returns the pressure [Pa] at CO2's triple point.
Definition CO2.hpp:137
static OPM_HOST_DEVICE Evaluation gasViscosity(const Params &params, Evaluation temperature, const Evaluation &pressure, bool extrapolate=false)
The dynamic viscosity [Pa s] of CO2.
Definition CO2.hpp:249
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:39
static constexpr Scalar R
The ideal gas constant .
Definition IdealGas.hpp:42
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30