Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PorousFlowBrineCO2.h"
11 : #include "BrineFluidProperties.h"
12 : #include "SinglePhaseFluidProperties.h"
13 : #include "MathUtils.h"
14 : #include "Conversion.h"
15 :
16 : registerMooseObject("PorousFlowApp", PorousFlowBrineCO2);
17 :
18 : InputParameters
19 938 : PorousFlowBrineCO2::validParams()
20 : {
21 938 : InputParameters params = PorousFlowFluidStateMultiComponentBase::validParams();
22 1876 : params.addRequiredParam<UserObjectName>("brine_fp", "The name of the user object for brine");
23 1876 : params.addRequiredParam<UserObjectName>("co2_fp", "The name of the user object for CO2");
24 1876 : params.addParam<unsigned int>("salt_component", 2, "The component number of salt");
25 938 : params.addClassDescription("Fluid state class for brine and CO2");
26 938 : return params;
27 0 : }
28 :
29 469 : PorousFlowBrineCO2::PorousFlowBrineCO2(const InputParameters & parameters)
30 : : PorousFlowFluidStateMultiComponentBase(parameters),
31 469 : _salt_component(getParam<unsigned int>("salt_component")),
32 469 : _brine_fp(getUserObject<BrineFluidProperties>("brine_fp")),
33 469 : _co2_fp(getUserObject<SinglePhaseFluidProperties>("co2_fp")),
34 469 : _Mh2o(_brine_fp.molarMassH2O()),
35 469 : _invMh2o(1.0 / _Mh2o),
36 469 : _Mco2(_co2_fp.molarMass()),
37 469 : _Mnacl(_brine_fp.molarMassNaCl()),
38 469 : _Rbar(_R * 10.0),
39 469 : _Tlower(372.15),
40 469 : _Tupper(382.15),
41 469 : _Zmin(1.0e-4),
42 938 : _co2_henry(_co2_fp.henryCoefficients())
43 : {
44 : // Check that the correct FluidProperties UserObjects have been provided
45 938 : if (_co2_fp.fluidName() != "co2")
46 0 : paramError("co2_fp", "A valid CO2 FluidProperties UserObject must be provided");
47 :
48 938 : if (_brine_fp.fluidName() != "brine")
49 0 : paramError("brine_fp", "A valid Brine FluidProperties UserObject must be provided");
50 :
51 : // Set the number of phases and components, and their indexes
52 469 : _num_phases = 2;
53 469 : _num_components = 3;
54 469 : _gas_phase_number = 1 - _aqueous_phase_number;
55 469 : _gas_fluid_component = 3 - _aqueous_fluid_component - _salt_component;
56 :
57 : // Check that _aqueous_phase_number is <= total number of phases
58 469 : if (_aqueous_phase_number >= _num_phases)
59 0 : paramError("liquid_phase_number",
60 : "This value is larger than the possible number of phases ",
61 : _num_phases);
62 :
63 : // Check that _aqueous_fluid_component is <= total number of fluid components
64 469 : if (_aqueous_fluid_component >= _num_components)
65 0 : paramError("liquid_fluid_component",
66 : "This value is larger than the possible number of fluid components",
67 : _num_components);
68 :
69 : // Check that the salt component index is not identical to the liquid fluid component
70 469 : if (_salt_component == _aqueous_fluid_component)
71 0 : paramError(
72 : "salt_component",
73 : "The value provided must be different from the value entered in liquid_fluid_component");
74 :
75 : // Check that _salt_component is <= total number of fluid components
76 469 : if (_salt_component >= _num_components)
77 0 : paramError("salt_component",
78 : "The value provided is larger than the possible number of fluid components",
79 : _num_components);
80 :
81 469 : _empty_fsp = FluidStateProperties(_num_components);
82 469 : }
83 :
84 : std::string
85 1 : PorousFlowBrineCO2::fluidStateName() const
86 : {
87 1 : return "brine-co2";
88 : }
89 :
90 : void
91 376495 : PorousFlowBrineCO2::thermophysicalProperties(Real pressure,
92 : Real temperature,
93 : Real Xnacl,
94 : Real Z,
95 : unsigned int qp,
96 : std::vector<FluidStateProperties> & fsp) const
97 : {
98 : // Make AD versions of primary variables then call AD thermophysicalProperties()
99 376495 : ADReal p = pressure;
100 376495 : Moose::derivInsert(p.derivatives(), _pidx, 1.0);
101 376495 : ADReal T = temperature;
102 376495 : Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
103 376495 : ADReal Zco2 = Z;
104 376495 : Moose::derivInsert(Zco2.derivatives(), _Zidx, 1.0);
105 376495 : ADReal X = Xnacl;
106 376495 : Moose::derivInsert(X.derivatives(), _Xidx, 1.0);
107 :
108 376495 : thermophysicalProperties(p, T, X, Zco2, qp, fsp);
109 376372 : }
110 :
111 : void
112 636106 : PorousFlowBrineCO2::thermophysicalProperties(const ADReal & pressure,
113 : const ADReal & temperature,
114 : const ADReal & Xnacl,
115 : const ADReal & Z,
116 : unsigned int qp,
117 : std::vector<FluidStateProperties> & fsp) const
118 : {
119 636106 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
120 636106 : FluidStateProperties & gas = fsp[_gas_phase_number];
121 :
122 : // Check whether the input temperature is within the region of validity
123 636106 : checkVariables(pressure.value(), temperature.value());
124 :
125 : // Clear all of the FluidStateProperties data
126 636106 : clearFluidStateProperties(fsp);
127 :
128 : FluidStatePhaseEnum phase_state;
129 636106 : massFractions(pressure, temperature, Xnacl, Z, phase_state, fsp);
130 :
131 636106 : switch (phase_state)
132 : {
133 8120 : case FluidStatePhaseEnum::GAS:
134 : {
135 : // Set the gas saturations
136 8120 : gas.saturation = 1.0;
137 :
138 : // Calculate gas properties
139 8120 : gasProperties(pressure, temperature, fsp);
140 :
141 : break;
142 : }
143 :
144 593207 : case FluidStatePhaseEnum::LIQUID:
145 : {
146 : // Calculate the liquid properties
147 593207 : const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0, qp);
148 593207 : liquidProperties(liquid_pressure, temperature, Xnacl, fsp);
149 :
150 : break;
151 : }
152 :
153 34779 : case FluidStatePhaseEnum::TWOPHASE:
154 : {
155 : // Calculate the gas and liquid properties in the two phase region
156 34779 : twoPhaseProperties(pressure, temperature, Xnacl, Z, qp, fsp);
157 :
158 : break;
159 : }
160 : }
161 :
162 : // Liquid saturations can now be set
163 1271966 : liquid.saturation = 1.0 - gas.saturation;
164 :
165 : // Save pressures to FluidStateProperties object
166 635983 : gas.pressure = pressure;
167 1271966 : liquid.pressure = pressure - _pc.capillaryPressure(liquid.saturation, qp);
168 635983 : }
169 :
170 : void
171 636135 : PorousFlowBrineCO2::massFractions(const ADReal & pressure,
172 : const ADReal & temperature,
173 : const ADReal & Xnacl,
174 : const ADReal & Z,
175 : FluidStatePhaseEnum & phase_state,
176 : std::vector<FluidStateProperties> & fsp) const
177 : {
178 636135 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
179 636135 : FluidStateProperties & gas = fsp[_gas_phase_number];
180 :
181 636135 : ADReal Xco2 = 0.0;
182 636135 : ADReal Yh2o = 0.0;
183 636135 : ADReal Yco2 = 0.0;
184 :
185 : // If the amount of CO2 is less than the smallest solubility, then all CO2 will
186 : // be dissolved, and the equilibrium mass fractions do not need to be computed
187 636135 : if (Z.value() < _Zmin)
188 501912 : phase_state = FluidStatePhaseEnum::LIQUID;
189 :
190 : else
191 : {
192 : // Equilibrium mass fraction of CO2 in liquid and H2O in gas phases
193 134223 : equilibriumMassFractions(pressure, temperature, Xnacl, Xco2, Yh2o);
194 :
195 268446 : Yco2 = 1.0 - Yh2o;
196 :
197 : // Determine which phases are present based on the value of z
198 134223 : phaseState(Z.value(), Xco2.value(), Yco2.value(), phase_state);
199 : }
200 :
201 : // The equilibrium mass fractions calculated above are only correct in the two phase
202 : // state. If only liquid or gas phases are present, the mass fractions are given by
203 : // the total mass fraction z
204 636135 : ADReal Xh2o = 0.0;
205 :
206 636135 : switch (phase_state)
207 : {
208 593211 : case FluidStatePhaseEnum::LIQUID:
209 : {
210 593211 : Xco2 = Z;
211 593211 : Yco2 = 0.0;
212 1186422 : Xh2o = 1.0 - Z;
213 593211 : Yh2o = 0.0;
214 593211 : break;
215 : }
216 :
217 8124 : case FluidStatePhaseEnum::GAS:
218 : {
219 8124 : Xco2 = 0.0;
220 8124 : Yco2 = Z;
221 16248 : Yh2o = 1.0 - Z;
222 8124 : break;
223 : }
224 :
225 34800 : case FluidStatePhaseEnum::TWOPHASE:
226 : {
227 : // Keep equilibrium mass fractions
228 69600 : Xh2o = 1.0 - Xco2;
229 34800 : break;
230 : }
231 : }
232 :
233 : // Save the mass fractions in the FluidStateProperties object
234 636135 : liquid.mass_fraction[_aqueous_fluid_component] = Xh2o;
235 636135 : liquid.mass_fraction[_gas_fluid_component] = Xco2;
236 636135 : liquid.mass_fraction[_salt_component] = Xnacl;
237 636135 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
238 636135 : gas.mass_fraction[_gas_fluid_component] = Yco2;
239 636135 : }
240 :
241 : void
242 42944 : PorousFlowBrineCO2::gasProperties(const ADReal & pressure,
243 : const ADReal & temperature,
244 : std::vector<FluidStateProperties> & fsp) const
245 : {
246 42944 : FluidStateProperties & gas = fsp[_gas_phase_number];
247 :
248 : // Gas density, viscosity and enthalpy are approximated with pure CO2 - no correction due
249 : // to the small amount of water vapor is made
250 : ADReal co2_density, co2_viscosity;
251 42944 : _co2_fp.rho_mu_from_p_T(pressure, temperature, co2_density, co2_viscosity);
252 :
253 42944 : ADReal co2_enthalpy = _co2_fp.h_from_p_T(pressure, temperature);
254 :
255 : // Save the values to the FluidStateProperties object. Note that derivatives wrt z are 0
256 42944 : gas.density = co2_density;
257 42944 : gas.viscosity = co2_viscosity;
258 42944 : gas.enthalpy = co2_enthalpy;
259 :
260 : mooseAssert(gas.density.value() > 0.0, "Gas density must be greater than zero");
261 42944 : gas.internal_energy = gas.enthalpy - pressure / gas.density;
262 42944 : }
263 :
264 : void
265 628042 : PorousFlowBrineCO2::liquidProperties(const ADReal & pressure,
266 : const ADReal & temperature,
267 : const ADReal & Xnacl,
268 : std::vector<FluidStateProperties> & fsp) const
269 : {
270 628042 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
271 :
272 : // The liquid density includes the density increase due to dissolved CO2
273 628042 : const ADReal brine_density = _brine_fp.rho_from_p_T_X(pressure, temperature, Xnacl);
274 :
275 : // Mass fraction of CO2 in liquid phase
276 627919 : const ADReal Xco2 = liquid.mass_fraction[_gas_fluid_component];
277 :
278 : // The liquid density
279 627919 : const ADReal co2_partial_density = partialDensityCO2(temperature);
280 :
281 1883757 : const ADReal liquid_density = 1.0 / (Xco2 / co2_partial_density + (1.0 - Xco2) / brine_density);
282 :
283 : // Assume that liquid viscosity is just the brine viscosity
284 627919 : const ADReal liquid_viscosity = _brine_fp.mu_from_p_T_X(pressure, temperature, Xnacl);
285 :
286 : // Liquid enthalpy (including contribution due to the enthalpy of dissolution)
287 627919 : const ADReal brine_enthalpy = _brine_fp.h_from_p_T_X(pressure, temperature, Xnacl);
288 :
289 : // Enthalpy of CO2
290 627919 : const ADReal co2_enthalpy = _co2_fp.h_from_p_T(pressure, temperature);
291 :
292 : // Enthalpy of dissolution
293 627919 : const ADReal hdis = enthalpyOfDissolution(temperature);
294 :
295 627919 : const ADReal liquid_enthalpy = (1.0 - Xco2) * brine_enthalpy + Xco2 * (co2_enthalpy + hdis);
296 :
297 : // Save the values to the FluidStateProperties object
298 627919 : liquid.density = liquid_density;
299 627919 : liquid.viscosity = liquid_viscosity;
300 627919 : liquid.enthalpy = liquid_enthalpy;
301 :
302 : mooseAssert(liquid.density.value() > 0.0, "Liquid density must be greater than zero");
303 627919 : liquid.internal_energy = liquid.enthalpy - pressure / liquid.density;
304 627919 : }
305 :
306 : ADReal
307 34790 : PorousFlowBrineCO2::saturation(const ADReal & pressure,
308 : const ADReal & temperature,
309 : const ADReal & Xnacl,
310 : const ADReal & Z,
311 : std::vector<FluidStateProperties> & fsp) const
312 : {
313 34790 : auto & gas = fsp[_gas_phase_number];
314 34790 : auto & liquid = fsp[_aqueous_fluid_component];
315 :
316 : // Approximate liquid density as saturation isn't known yet, by using the gas
317 : // pressure rather than the liquid pressure. This does result in a small error
318 : // in the calculated saturation, but this is below the error associated with
319 : // the correlations. A more accurate saturation could be found iteraviely,
320 : // at the cost of increased computational expense
321 :
322 : // Gas density
323 34790 : const ADReal gas_density = _co2_fp.rho_from_p_T(pressure, temperature);
324 :
325 : // Approximate liquid density as saturation isn't known yet
326 34790 : const ADReal brine_density = _brine_fp.rho_from_p_T_X(pressure, temperature, Xnacl);
327 :
328 : // Mass fraction of CO2 in liquid phase
329 34790 : const ADReal Xco2 = liquid.mass_fraction[_gas_fluid_component];
330 :
331 : // The liquid density
332 34790 : const ADReal co2_partial_density = partialDensityCO2(temperature);
333 :
334 104370 : const ADReal liquid_density = 1.0 / (Xco2 / co2_partial_density + (1.0 - Xco2) / brine_density);
335 :
336 34790 : const ADReal Yco2 = gas.mass_fraction[_gas_fluid_component];
337 :
338 : // Set mass equilibrium constants used in the calculation of vapor mass fraction
339 : const ADReal K0 = Yco2 / Xco2;
340 69580 : const ADReal K1 = (1.0 - Yco2) / (1.0 - Xco2);
341 34790 : const ADReal vapor_mass_fraction = vaporMassFraction(Z, K0, K1);
342 :
343 : // The gas saturation in the two phase case
344 34790 : const ADReal saturation = vapor_mass_fraction * liquid_density /
345 0 : (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
346 :
347 34790 : return saturation;
348 : }
349 :
350 : void
351 34779 : PorousFlowBrineCO2::twoPhaseProperties(const ADReal & pressure,
352 : const ADReal & temperature,
353 : const ADReal & Xnacl,
354 : const ADReal & Z,
355 : unsigned int qp,
356 : std::vector<FluidStateProperties> & fsp) const
357 : {
358 34779 : auto & gas = fsp[_gas_phase_number];
359 :
360 : // Calculate all of the gas phase properties, as these don't depend on saturation
361 34779 : gasProperties(pressure, temperature, fsp);
362 :
363 : // The gas saturation in the two phase case
364 34779 : gas.saturation = saturation(pressure, temperature, Xnacl, Z, fsp);
365 :
366 : // The liquid pressure and properties can now be calculated
367 69558 : const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0 - gas.saturation, qp);
368 34779 : liquidProperties(liquid_pressure, temperature, Xnacl, fsp);
369 34779 : }
370 :
371 : void
372 134269 : PorousFlowBrineCO2::equilibriumMassFractions(const ADReal & pressure,
373 : const ADReal & temperature,
374 : const ADReal & Xnacl,
375 : ADReal & Xco2,
376 : ADReal & Yh2o) const
377 : {
378 : // Mole fractions at equilibrium
379 : ADReal xco2, yh2o;
380 134269 : equilibriumMoleFractions(pressure, temperature, Xnacl, xco2, yh2o);
381 :
382 : // The mass fraction of H2O in gas (assume no salt in gas phase) and derivatives
383 : // wrt p, T, and X
384 402807 : Yh2o = yh2o * _Mh2o / (yh2o * _Mh2o + (1.0 - yh2o) * _Mco2);
385 :
386 : // NaCl molality (mol/kg)
387 134269 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
388 :
389 : // The molality of CO2 in 1kg of H2O
390 402807 : const ADReal mco2 = xco2 * (2.0 * mnacl + _invMh2o) / (1.0 - xco2);
391 : // The mass fraction of CO2 in brine is then
392 268538 : const ADReal denominator = (1.0 + mnacl * _Mnacl + mco2 * _Mco2);
393 134269 : Xco2 = mco2 * _Mco2 / denominator;
394 134269 : }
395 :
396 : void
397 133865 : PorousFlowBrineCO2::fugacityCoefficientsLowTemp(const ADReal & pressure,
398 : const ADReal & temperature,
399 : const ADReal & co2_density,
400 : ADReal & fco2,
401 : ADReal & fh2o) const
402 : {
403 : using std::log, std::exp, std::pow;
404 :
405 133865 : if (temperature.value() > 373.15)
406 0 : mooseError(name(),
407 : ": fugacityCoefficientsLowTemp() is not valid for T > 373.15K. Use "
408 : "fugacityCoefficientsHighTemp() instead");
409 :
410 : // Need pressure in bar
411 133865 : const ADReal pbar = pressure * 1.0e-5;
412 :
413 : // Molar volume in cm^3/mol
414 133865 : const ADReal V = _Mco2 / co2_density * 1.0e6;
415 :
416 : // Redlich-Kwong parameters
417 401595 : const ADReal aCO2 = 7.54e7 - 4.13e4 * temperature;
418 133865 : const Real bCO2 = 27.8;
419 133865 : const Real aCO2H2O = 7.89e7;
420 133865 : const Real bH2O = 18.18;
421 :
422 133865 : const ADReal t15 = pow(temperature, 1.5);
423 :
424 : // The fugacity coefficients for H2O and CO2
425 535460 : auto lnPhi = [V, aCO2, bCO2, t15, this](ADReal a, ADReal b)
426 : {
427 267730 : return log(V / (V - bCO2)) + b / (V - bCO2) -
428 535460 : 2.0 * a / (_Rbar * t15 * bCO2) * log((V + bCO2) / V) +
429 1338650 : aCO2 * b / (_Rbar * t15 * bCO2 * bCO2) * (log((V + bCO2) / V) - bCO2 / (V + bCO2));
430 133865 : };
431 :
432 401595 : const ADReal lnPhiH2O = lnPhi(aCO2H2O, bH2O) - log(pbar * V / (_Rbar * temperature));
433 133865 : const ADReal lnPhiCO2 = lnPhi(aCO2, bCO2) - log(pbar * V / (_Rbar * temperature));
434 :
435 133865 : fh2o = exp(lnPhiH2O);
436 133865 : fco2 = exp(lnPhiCO2);
437 133865 : }
438 :
439 : void
440 10402 : PorousFlowBrineCO2::fugacityCoefficientsHighTemp(const ADReal & pressure,
441 : const ADReal & temperature,
442 : const ADReal & co2_density,
443 : const ADReal & xco2,
444 : const ADReal & yh2o,
445 : ADReal & fco2,
446 : ADReal & fh2o) const
447 : {
448 10402 : if (temperature.value() <= 373.15)
449 0 : mooseError(name(),
450 : ": fugacityCoefficientsHighTemp() is not valid for T <= 373.15K. Use "
451 : "fugacityCoefficientsLowTemp() instead");
452 :
453 10402 : fh2o = fugacityCoefficientH2OHighTemp(pressure, temperature, co2_density, xco2, yh2o);
454 10402 : fco2 = fugacityCoefficientCO2HighTemp(pressure, temperature, co2_density, xco2, yh2o);
455 10402 : }
456 :
457 : ADReal
458 10403 : PorousFlowBrineCO2::fugacityCoefficientH2OHighTemp(const ADReal & pressure,
459 : const ADReal & temperature,
460 : const ADReal & co2_density,
461 : const ADReal & xco2,
462 : const ADReal & yh2o) const
463 : {
464 : using std::sqrt, std::pow, std::log, std::exp;
465 :
466 : // Need pressure in bar
467 10403 : const ADReal pbar = pressure * 1.0e-5;
468 : // Molar volume in cm^3/mol
469 10403 : const ADReal V = _Mco2 / co2_density * 1.0e6;
470 :
471 : // Redlich-Kwong parameters
472 20806 : const ADReal yco2 = 1.0 - yh2o;
473 10403 : const ADReal xh2o = 1.0 - xco2;
474 :
475 31209 : const ADReal aCO2 = 8.008e7 - 4.984e4 * temperature;
476 31209 : const ADReal aH2O = 1.337e8 - 1.4e4 * temperature;
477 10403 : const Real bCO2 = 28.25;
478 10403 : const Real bH2O = 15.7;
479 31209 : const ADReal KH2OCO2 = 1.427e-2 - 4.037e-4 * temperature;
480 31209 : const ADReal KCO2H2O = 0.4228 - 7.422e-4 * temperature;
481 10403 : const ADReal kH2OCO2 = KH2OCO2 * yh2o + KCO2H2O * yco2;
482 10403 : const ADReal kCO2H2O = KCO2H2O * yh2o + KH2OCO2 * yco2;
483 :
484 20806 : const ADReal aH2OCO2 = sqrt(aCO2 * aH2O) * (1.0 - kH2OCO2);
485 31209 : const ADReal aCO2H2O = sqrt(aCO2 * aH2O) * (1.0 - kCO2H2O);
486 :
487 10403 : const ADReal amix = yh2o * yh2o * aH2O + yh2o * yco2 * (aH2OCO2 + aCO2H2O) + yco2 * yco2 * aCO2;
488 10403 : const ADReal bmix = yh2o * bH2O + yco2 * bCO2;
489 :
490 10403 : const ADReal t15 = pow(temperature, 1.5);
491 :
492 0 : ADReal lnPhiH2O = bH2O / bmix * (pbar * V / (_Rbar * temperature) - 1.0) -
493 31209 : log(pbar * (V - bmix) / (_Rbar * temperature));
494 10403 : ADReal term3 = (2.0 * yh2o * aH2O + yco2 * (aH2OCO2 + aCO2H2O) -
495 10403 : yh2o * yco2 * sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O) * (yh2o - yco2) +
496 20806 : xh2o * xco2 * sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O)) /
497 : amix;
498 10403 : term3 -= bH2O / bmix;
499 20806 : term3 *= amix / (bmix * _Rbar * t15) * log(V / (V + bmix));
500 10403 : lnPhiH2O += term3;
501 :
502 10403 : return exp(lnPhiH2O);
503 : }
504 :
505 : ADReal
506 10403 : PorousFlowBrineCO2::fugacityCoefficientCO2HighTemp(const ADReal & pressure,
507 : const ADReal & temperature,
508 : const ADReal & co2_density,
509 : const ADReal & xco2,
510 : const ADReal & yh2o) const
511 : {
512 : using std::sqrt, std::pow, std::log, std::exp;
513 :
514 : // Need pressure in bar
515 10403 : const ADReal pbar = pressure * 1.0e-5;
516 : // Molar volume in cm^3/mol
517 10403 : const ADReal V = _Mco2 / co2_density * 1.0e6;
518 :
519 : // Redlich-Kwong parameters
520 10403 : const ADReal yco2 = 1.0 - yh2o;
521 10403 : const ADReal xh2o = 1.0 - xco2;
522 :
523 31209 : const ADReal aCO2 = 8.008e7 - 4.984e4 * temperature;
524 31209 : const ADReal aH2O = 1.337e8 - 1.4e4 * temperature;
525 10403 : const Real bCO2 = 28.25;
526 10403 : const Real bH2O = 15.7;
527 31209 : const ADReal KH2OCO2 = 1.427e-2 - 4.037e-4 * temperature;
528 31209 : const ADReal KCO2H2O = 0.4228 - 7.422e-4 * temperature;
529 10403 : const ADReal kH2OCO2 = KH2OCO2 * yh2o + KCO2H2O * yco2;
530 10403 : const ADReal kCO2H2O = KCO2H2O * yh2o + KH2OCO2 * yco2;
531 :
532 20806 : const ADReal aH2OCO2 = sqrt(aCO2 * aH2O) * (1.0 - kH2OCO2);
533 31209 : const ADReal aCO2H2O = sqrt(aCO2 * aH2O) * (1.0 - kCO2H2O);
534 :
535 10403 : const ADReal amix = yh2o * yh2o * aH2O + yh2o * yco2 * (aH2OCO2 + aCO2H2O) + yco2 * yco2 * aCO2;
536 10403 : const ADReal bmix = yh2o * bH2O + yco2 * bCO2;
537 :
538 10403 : const ADReal t15 = pow(temperature, 1.5);
539 :
540 0 : ADReal lnPhiCO2 = bCO2 / bmix * (pbar * V / (_Rbar * temperature) - 1.0) -
541 31209 : log(pbar * (V - bmix) / (_Rbar * temperature));
542 :
543 10403 : ADReal term3 = (2.0 * yco2 * aCO2 + yh2o * (aH2OCO2 + aCO2H2O) -
544 10403 : yh2o * yco2 * sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O) * (yh2o - yco2) +
545 20806 : xh2o * xco2 * sqrt(aH2O * aCO2) * (kCO2H2O - kH2OCO2)) /
546 : amix;
547 :
548 20806 : lnPhiCO2 += (term3 - bCO2 / bmix) * amix / (bmix * _Rbar * t15) * log(V / (V + bmix));
549 :
550 10403 : return exp(lnPhiCO2);
551 : }
552 :
553 : ADReal
554 10403 : PorousFlowBrineCO2::activityCoefficientH2O(const ADReal & temperature, const ADReal & xco2) const
555 : {
556 : using std::exp;
557 :
558 10403 : if (temperature.value() <= 373.15)
559 1 : return 1.0;
560 : else
561 : {
562 : const ADReal Tref = temperature - 373.15;
563 10402 : const ADReal xh2o = 1.0 - xco2;
564 20804 : const ADReal Am = -3.084e-2 * Tref + 1.927e-5 * Tref * Tref;
565 :
566 20804 : return exp((Am - 2.0 * Am * xh2o) * xco2 * xco2);
567 : }
568 : }
569 :
570 : ADReal
571 10403 : PorousFlowBrineCO2::activityCoefficientCO2(const ADReal & temperature, const ADReal & xco2) const
572 : {
573 : using std::exp;
574 :
575 10403 : if (temperature.value() <= 373.15)
576 1 : return 1.0;
577 : else
578 : {
579 : const ADReal Tref = temperature - 373.15;
580 10402 : const ADReal xh2o = 1.0 - xco2;
581 20804 : const ADReal Am = -3.084e-2 * Tref + 1.927e-5 * Tref * Tref;
582 :
583 20804 : return exp(2.0 * Am * xco2 * xh2o * xh2o);
584 : }
585 : }
586 :
587 : ADReal
588 133869 : PorousFlowBrineCO2::activityCoefficient(const ADReal & pressure,
589 : const ADReal & temperature,
590 : const ADReal & Xnacl) const
591 : {
592 : using std::log, std::exp;
593 :
594 : // Need pressure in bar
595 133869 : const ADReal pbar = pressure * 1.0e-5;
596 : // Need NaCl molality (mol/kg)
597 133869 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
598 :
599 401607 : const ADReal lambda = -0.411370585 + 6.07632013e-4 * temperature + 97.5347708 / temperature -
600 133869 : 0.0237622469 * pbar / temperature +
601 267738 : 0.0170656236 * pbar / (630.0 - temperature) +
602 133869 : 1.41335834e-5 * temperature * log(pbar);
603 :
604 267738 : const ADReal xi = 3.36389723e-4 - 1.9829898e-5 * temperature +
605 133869 : 2.12220830e-3 * pbar / temperature -
606 401607 : 5.24873303e-3 * pbar / (630.0 - temperature);
607 :
608 267738 : return exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
609 : }
610 :
611 : ADReal
612 10405 : PorousFlowBrineCO2::activityCoefficientHighTemp(const ADReal & temperature,
613 : const ADReal & Xnacl) const
614 : {
615 : using std::exp;
616 :
617 : // Need NaCl molality (mol/kg)
618 20810 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
619 :
620 : const ADReal T2 = temperature * temperature;
621 : const ADReal T3 = temperature * T2;
622 :
623 31215 : const ADReal lambda = 2.217e-4 * temperature + 1.074 / temperature + 2648.0 / T2;
624 41620 : const ADReal xi = 1.3e-5 * temperature - 20.12 / temperature + 5259.0 / T2;
625 :
626 41620 : return (1.0 - mnacl / _invMh2o) * exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
627 : }
628 :
629 : ADReal
630 144269 : PorousFlowBrineCO2::equilibriumConstantH2O(const ADReal & temperature) const
631 : {
632 : using std::pow;
633 :
634 : // Uses temperature in Celsius
635 : const ADReal Tc = temperature - _T_c2k;
636 : const ADReal Tc2 = Tc * Tc;
637 : const ADReal Tc3 = Tc2 * Tc;
638 : const ADReal Tc4 = Tc3 * Tc;
639 :
640 : ADReal logK0H2O;
641 :
642 144269 : if (Tc <= 99.0)
643 669330 : logK0H2O = -2.209 + 3.097e-2 * Tc - 1.098e-4 * Tc2 + 2.048e-7 * Tc3;
644 :
645 10403 : else if (Tc > 99.0 && Tc < 109.0)
646 : {
647 0 : const ADReal Tint = (Tc - 99.0) / 10.0;
648 : const ADReal Tint2 = Tint * Tint;
649 0 : logK0H2O = -0.0204026 + 0.0152513 * Tint + 0.417565 * Tint2 - 0.278636 * Tint * Tint2;
650 : }
651 :
652 : else // 109 <= Tc <= 300
653 62418 : logK0H2O = -2.1077 + 2.8127e-2 * Tc - 8.4298e-5 * Tc2 + 1.4969e-7 * Tc3 - 1.1812e-10 * Tc4;
654 :
655 144269 : return pow(10.0, logK0H2O);
656 : }
657 :
658 : ADReal
659 144269 : PorousFlowBrineCO2::equilibriumConstantCO2(const ADReal & temperature) const
660 : {
661 : using std::pow;
662 :
663 : // Uses temperature in Celsius
664 : const ADReal Tc = temperature - _T_c2k;
665 : const ADReal Tc2 = Tc * Tc;
666 : const ADReal Tc3 = Tc2 * Tc;
667 :
668 : ADReal logK0CO2;
669 :
670 144269 : if (Tc <= 99.0)
671 535464 : logK0CO2 = 1.189 + 1.304e-2 * Tc - 5.446e-5 * Tc2;
672 :
673 10403 : else if (Tc > 99.0 && Tc < 109.0)
674 : {
675 0 : const ADReal Tint = (Tc - 99.0) / 10.0;
676 : const ADReal Tint2 = Tint * Tint;
677 0 : logK0CO2 = 1.9462 + 2.25692e-2 * Tint - 9.49577e-3 * Tint2 - 6.77721e-3 * Tint * Tint2;
678 : }
679 :
680 : else // 109 <= Tc <= 300
681 52015 : logK0CO2 = 1.668 + 3.992e-3 * Tc - 1.156e-5 * Tc2 + 1.593e-9 * Tc3;
682 :
683 144269 : return pow(10.0, logK0CO2);
684 : }
685 :
686 : void
687 134276 : PorousFlowBrineCO2::equilibriumMoleFractions(const ADReal & pressure,
688 : const ADReal & temperature,
689 : const ADReal & Xnacl,
690 : ADReal & xco2,
691 : ADReal & yh2o) const
692 : {
693 134276 : if (temperature.value() <= _Tlower)
694 : {
695 133862 : equilibriumMoleFractionsLowTemp(pressure, temperature, Xnacl, xco2, yh2o);
696 : }
697 414 : else if (temperature.value() > _Tlower && temperature.value() < _Tupper)
698 : {
699 : // Cubic polynomial in this regime
700 2 : const Real Tint = (temperature.value() - _Tlower) / 10.0;
701 :
702 : // Equilibrium mole fractions and derivatives at the lower temperature
703 2 : ADReal Tlower = _Tlower;
704 2 : Moose::derivInsert(Tlower.derivatives(), _Tidx, 1.0);
705 :
706 : ADReal xco2_lower, yh2o_lower;
707 2 : equilibriumMoleFractionsLowTemp(pressure, Tlower, Xnacl, xco2_lower, yh2o_lower);
708 :
709 2 : const Real dxco2_dT_lower = xco2_lower.derivatives()[_Tidx];
710 2 : const Real dyh2o_dT_lower = yh2o_lower.derivatives()[_Tidx];
711 :
712 : // Equilibrium mole fractions and derivatives at the upper temperature
713 : Real xco2_upper, yh2o_upper;
714 2 : Real co2_density_upper = _co2_fp.rho_from_p_T(pressure.value(), _Tupper);
715 :
716 2 : solveEquilibriumMoleFractionHighTemp(
717 2 : pressure.value(), _Tupper, Xnacl.value(), co2_density_upper, xco2_upper, yh2o_upper);
718 :
719 : Real A, dA_dp, dA_dT, B, dB_dp, dB_dT, dB_dX;
720 2 : funcABHighTemp(pressure.value(),
721 2 : _Tupper,
722 : Xnacl.value(),
723 : co2_density_upper,
724 : xco2_upper,
725 : yh2o_upper,
726 : A,
727 : dA_dp,
728 : dA_dT,
729 : B,
730 : dB_dp,
731 : dB_dT,
732 : dB_dX);
733 :
734 2 : const Real dyh2o_dT_upper =
735 2 : ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
736 2 : const Real dxco2_dT_upper = dB_dT * (1.0 - yh2o_upper) - B * dyh2o_dT_upper;
737 :
738 : // The mole fractions in this regime are then found by interpolation
739 : Real xco2r, yh2or, dxco2_dT, dyh2o_dT;
740 2 : smoothCubicInterpolation(
741 : Tint, xco2_lower.value(), dxco2_dT_lower, xco2_upper, dxco2_dT_upper, xco2r, dxco2_dT);
742 2 : smoothCubicInterpolation(
743 : Tint, yh2o_lower.value(), dyh2o_dT_lower, yh2o_upper, dyh2o_dT_upper, yh2or, dyh2o_dT);
744 :
745 2 : xco2 = xco2r;
746 4 : Moose::derivInsert(xco2.derivatives(), _pidx, xco2_lower.derivatives()[_pidx]);
747 2 : Moose::derivInsert(xco2.derivatives(), _Tidx, dxco2_dT);
748 4 : Moose::derivInsert(xco2.derivatives(), _Xidx, xco2_lower.derivatives()[_Xidx]);
749 :
750 2 : yh2o = yh2or;
751 4 : Moose::derivInsert(yh2o.derivatives(), _pidx, yh2o_lower.derivatives()[_pidx]);
752 2 : Moose::derivInsert(yh2o.derivatives(), _Tidx, dyh2o_dT);
753 4 : Moose::derivInsert(yh2o.derivatives(), _Xidx, yh2o_lower.derivatives()[_Xidx]);
754 : }
755 : else
756 : {
757 : // CO2 density and derivatives wrt pressure and temperature
758 412 : const Real co2_density = _co2_fp.rho_from_p_T(pressure.value(), temperature.value());
759 :
760 : // Equilibrium mole fractions solved using iteration in this regime
761 : Real xco2r, yh2or;
762 412 : solveEquilibriumMoleFractionHighTemp(
763 : pressure.value(), temperature.value(), Xnacl.value(), co2_density, xco2r, yh2or);
764 :
765 : // Can use these in funcABHighTemp() to get derivatives analytically rather than by iteration
766 : Real A, dA_dp, dA_dT, B, dB_dp, dB_dT, dB_dX;
767 412 : funcABHighTemp(pressure.value(),
768 : temperature.value(),
769 : Xnacl.value(),
770 : co2_density,
771 : xco2r,
772 : yh2or,
773 : A,
774 : dA_dp,
775 : dA_dT,
776 : B,
777 : dB_dp,
778 : dB_dT,
779 : dB_dX);
780 :
781 412 : const Real dyh2o_dp =
782 412 : ((1.0 - B) * dA_dp + (A - 1.0) * A * dB_dp) / (1.0 - A * B) / (1.0 - A * B);
783 412 : const Real dxco2_dp = dB_dp * (1.0 - yh2or) - B * dyh2o_dp;
784 :
785 412 : const Real dyh2o_dT =
786 412 : ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
787 412 : const Real dxco2_dT = dB_dT * (1.0 - yh2or) - B * dyh2o_dT;
788 :
789 412 : const Real dyh2o_dX = ((A - 1.0) * A * dB_dX) / (1.0 - A * B) / (1.0 - A * B);
790 412 : const Real dxco2_dX = dB_dX * (1.0 - yh2or) - B * dyh2o_dX;
791 :
792 412 : xco2 = xco2r;
793 412 : Moose::derivInsert(xco2.derivatives(), _pidx, dxco2_dp);
794 412 : Moose::derivInsert(xco2.derivatives(), _Tidx, dxco2_dT);
795 412 : Moose::derivInsert(xco2.derivatives(), _Xidx, dxco2_dX);
796 :
797 412 : yh2o = yh2or;
798 412 : Moose::derivInsert(yh2o.derivatives(), _pidx, dyh2o_dp);
799 412 : Moose::derivInsert(yh2o.derivatives(), _Tidx, dyh2o_dT);
800 412 : Moose::derivInsert(yh2o.derivatives(), _Xidx, dyh2o_dX);
801 : }
802 134276 : }
803 :
804 : void
805 133864 : PorousFlowBrineCO2::equilibriumMoleFractionsLowTemp(const ADReal & pressure,
806 : const ADReal & temperature,
807 : const ADReal & Xnacl,
808 : ADReal & xco2,
809 : ADReal & yh2o) const
810 : {
811 133864 : if (temperature.value() > 373.15)
812 0 : mooseError(name(),
813 : ": equilibriumMoleFractionsLowTemp() is not valid for T > 373.15K. Use "
814 : "equilibriumMoleFractions() instead");
815 :
816 : // CO2 density and derivatives wrt pressure and temperature
817 133864 : const ADReal co2_density = _co2_fp.rho_from_p_T(pressure, temperature);
818 :
819 : // Assume infinite dilution (yh20 = 0 and xco2 = 0) in low temperature regime
820 : ADReal A, B;
821 133864 : funcABLowTemp(pressure, temperature, co2_density, A, B);
822 :
823 : // As the activity coefficient for CO2 in brine used in this regime isn't a 'true'
824 : // activity coefficient, we instead calculate the molality of CO2 in water, then
825 : // correct it for brine, and then calculate the mole fractions.
826 : // The mole fraction in pure water is
827 267728 : const ADReal yh2ow = (1.0 - B) / (1.0 / A - B);
828 133864 : const ADReal xco2w = B * (1.0 - yh2ow);
829 :
830 : // Molality of CO2 in pure water
831 267728 : const ADReal mco2w = xco2w * _invMh2o / (1.0 - xco2w);
832 : // Molality of CO2 in brine is then calculated using gamma
833 133864 : const ADReal gamma = activityCoefficient(pressure, temperature, Xnacl);
834 : const ADReal mco2 = mco2w / gamma;
835 :
836 : // Need NaCl molality (mol/kg)
837 133864 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
838 :
839 : // Mole fractions of CO2 and H2O in liquid and gas phases
840 267728 : const ADReal total_moles = 2.0 * mnacl + _invMh2o + mco2;
841 133864 : xco2 = mco2 / total_moles;
842 401592 : yh2o = A * (1.0 - xco2 - 2.0 * mnacl / total_moles);
843 133864 : }
844 :
845 : void
846 133864 : PorousFlowBrineCO2::funcABLowTemp(const ADReal & pressure,
847 : const ADReal & temperature,
848 : const ADReal & co2_density,
849 : ADReal & A,
850 : ADReal & B) const
851 : {
852 : using std::exp;
853 :
854 133864 : if (temperature.value() > 373.15)
855 0 : mooseError(name(),
856 : ": funcABLowTemp() is not valid for T > 373.15K. Use funcABHighTemp() instead");
857 :
858 : // Pressure in bar
859 133864 : const ADReal pbar = pressure * 1.0e-5;
860 :
861 : // Reference pressure and partial molar volumes
862 : const Real pref = 1.0;
863 133864 : const Real vCO2 = 32.6;
864 133864 : const Real vH2O = 18.1;
865 :
866 : const ADReal delta_pbar = pbar - pref;
867 133864 : const ADReal Rt = _Rbar * temperature;
868 :
869 : // Equilibrium constants
870 133864 : const ADReal K0H2O = equilibriumConstantH2O(temperature);
871 133864 : const ADReal K0CO2 = equilibriumConstantCO2(temperature);
872 :
873 : // Fugacity coefficients
874 : ADReal phiH2O, phiCO2;
875 133864 : fugacityCoefficientsLowTemp(pressure, temperature, co2_density, phiCO2, phiH2O);
876 :
877 267728 : A = K0H2O / (phiH2O * pbar) * exp(delta_pbar * vH2O / Rt);
878 401592 : B = phiCO2 * pbar / (_invMh2o * K0CO2) * exp(-delta_pbar * vCO2 / Rt);
879 133864 : }
880 :
881 : void
882 10401 : PorousFlowBrineCO2::funcABHighTemp(Real pressure,
883 : Real temperature,
884 : Real Xnacl,
885 : Real co2_density,
886 : Real xco2,
887 : Real yh2o,
888 : Real & A,
889 : Real & B) const
890 : {
891 : using std::exp;
892 :
893 10401 : if (temperature <= 373.15)
894 0 : mooseError(name(),
895 : ": funcABHighTemp() is not valid for T <= 373.15K. Use funcABLowTemp() instead");
896 :
897 : // Pressure in bar
898 10401 : const Real pbar = pressure * 1.0e-5;
899 : // Temperature in C
900 10401 : const Real Tc = temperature - _T_c2k;
901 :
902 : // Reference pressure and partial molar volumes
903 10401 : const Real pref = -1.9906e-1 + 2.0471e-3 * Tc + 1.0152e-4 * Tc * Tc - 1.4234e-6 * Tc * Tc * Tc +
904 10401 : 1.4168e-8 * Tc * Tc * Tc * Tc;
905 10401 : const Real vCO2 = 32.6 + 3.413e-2 * (Tc - 100.0);
906 10401 : const Real vH2O = 18.1 + 3.137e-2 * (Tc - 100.0);
907 :
908 10401 : const Real delta_pbar = pbar - pref;
909 10401 : const Real Rt = _Rbar * temperature;
910 :
911 : // Equilibrium constants
912 : // Use dummy ADReal temperature as derivatives aren't required
913 10401 : const ADReal T = temperature;
914 10401 : Real K0H2O = equilibriumConstantH2O(T).value();
915 10401 : Real K0CO2 = equilibriumConstantCO2(T).value();
916 :
917 : // Fugacity coefficients
918 : // Use dummy ADReal variables as derivatives aren't required
919 10401 : const ADReal p = pressure;
920 10401 : const ADReal rhoco2 = co2_density;
921 10401 : const ADReal x = xco2;
922 10401 : const ADReal y = yh2o;
923 :
924 : ADReal phiH2O, phiCO2;
925 10401 : fugacityCoefficientsHighTemp(p, T, rhoco2, x, y, phiCO2, phiH2O);
926 :
927 : // Activity coefficients
928 10401 : const Real gammaH2O = activityCoefficientH2O(T, x).value();
929 10401 : const Real gammaCO2 = activityCoefficientCO2(T, x).value();
930 :
931 : // Activity coefficient for CO2 in brine
932 : // Use dummy ADReal Xnacl as derivatives aren't required
933 10401 : const ADReal X = Xnacl;
934 10401 : const Real gamma = activityCoefficientHighTemp(T, X).value();
935 :
936 10401 : A = K0H2O * gammaH2O / (phiH2O.value() * pbar) * exp(delta_pbar * vH2O / Rt);
937 10401 : B = phiCO2.value() * pbar / (_invMh2o * K0CO2 * gamma * gammaCO2) * exp(-delta_pbar * vCO2 / Rt);
938 10401 : }
939 :
940 : void
941 414 : PorousFlowBrineCO2::funcABHighTemp(Real pressure,
942 : Real temperature,
943 : Real Xnacl,
944 : Real co2_density,
945 : Real xco2,
946 : Real yh2o,
947 : Real & A,
948 : Real & dA_dp,
949 : Real & dA_dT,
950 : Real & B,
951 : Real & dB_dp,
952 : Real & dB_dT,
953 : Real & dB_dX) const
954 : {
955 414 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, xco2, yh2o, A, B);
956 :
957 : // Use finite differences for derivatives in the high temperature regime
958 : const Real dp = 1.0e-2;
959 : const Real dT = 1.0e-6;
960 : const Real dX = 1.0e-8;
961 :
962 : Real A2, B2;
963 414 : funcABHighTemp(pressure + dp, temperature, Xnacl, co2_density, xco2, yh2o, A2, B2);
964 414 : dA_dp = (A2 - A) / dp;
965 414 : dB_dp = (B2 - B) / dp;
966 :
967 414 : funcABHighTemp(pressure, temperature + dT, Xnacl, co2_density, xco2, yh2o, A2, B2);
968 414 : dA_dT = (A2 - A) / dT;
969 414 : dB_dT = (B2 - B) / dT;
970 :
971 414 : funcABHighTemp(pressure, temperature, Xnacl + dX, co2_density, xco2, yh2o, A2, B2);
972 414 : dB_dX = (B2 - B) / dX;
973 414 : }
974 :
975 : void
976 418 : PorousFlowBrineCO2::solveEquilibriumMoleFractionHighTemp(
977 : Real pressure, Real temperature, Real Xnacl, Real co2_density, Real & xco2, Real & yh2o) const
978 : {
979 : // Initial guess for yh2o and xco2 (from Spycher and Pruess (2010))
980 418 : Real y = _brine_fp.vaporPressure(temperature, 0.0) / pressure;
981 : Real x = 0.009;
982 :
983 : // Need salt mass fraction in molality
984 418 : const Real mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
985 :
986 : // If y > 1, then just use y = 1, x = 0 (only a gas phase)
987 418 : if (y >= 1.0)
988 : {
989 : y = 1.0;
990 : x = 0.0;
991 : }
992 : else
993 : {
994 : // Residual function for Newton-Raphson
995 8745 : auto fy = [mnacl, this](Real y, Real A, Real B) {
996 : return y -
997 8745 : (1.0 - B) * _invMh2o / ((1.0 / A - B) * (2.0 * mnacl + _invMh2o) + 2.0 * mnacl * B);
998 417 : };
999 :
1000 : // Derivative of fy wrt y
1001 4164 : auto dfy = [mnacl, this](Real A, Real B, Real dA, Real dB)
1002 : {
1003 4164 : const Real denominator = (1.0 / A - B) * (2.0 * mnacl + _invMh2o) + 2.0 * mnacl * B;
1004 4164 : return 1.0 + _invMh2o * dB / denominator +
1005 4164 : (1.0 - B) * _invMh2o *
1006 4164 : (2.0 * mnacl * dB - (2.0 * mnacl + _invMh2o) * (dB + dA / A / A)) / denominator /
1007 4164 : denominator;
1008 417 : };
1009 :
1010 : Real A, B;
1011 : Real dA, dB;
1012 : const Real dy = 1.0e-8;
1013 :
1014 : // Solve for yh2o using Newton-Raphson method
1015 : unsigned int iter = 0;
1016 : const Real tol = 1.0e-12;
1017 : const unsigned int max_its = 10;
1018 417 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y, A, B);
1019 :
1020 4581 : while (std::abs(fy(y, A, B)) > tol)
1021 : {
1022 4164 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y, A, B);
1023 : // Finite difference derivatives of A and B wrt y
1024 4164 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y + dy, dA, dB);
1025 4164 : dA = (dA - A) / dy;
1026 4164 : dB = (dB - B) / dy;
1027 :
1028 4164 : y = y - fy(y, A, B) / dfy(A, B, dA, dB);
1029 :
1030 4164 : x = B * (1.0 - y);
1031 :
1032 : // Break if not converged and just use the value
1033 : if (iter > max_its)
1034 : break;
1035 : }
1036 : }
1037 :
1038 418 : yh2o = y;
1039 418 : xco2 = x;
1040 418 : }
1041 :
1042 : ADReal
1043 662712 : PorousFlowBrineCO2::partialDensityCO2(const ADReal & temperature) const
1044 : {
1045 : // This correlation uses temperature in C
1046 : const ADReal Tc = temperature - _T_c2k;
1047 : // The parial molar volume
1048 2650848 : const ADReal V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
1049 :
1050 1325424 : return 1.0e6 * _Mco2 / V;
1051 : }
1052 :
1053 : Real
1054 37 : PorousFlowBrineCO2::totalMassFraction(
1055 : Real pressure, Real temperature, Real Xnacl, Real saturation, unsigned int qp) const
1056 : {
1057 : // Check whether the input pressure and temperature are within the region of validity
1058 37 : checkVariables(pressure, temperature);
1059 :
1060 : // As we do not require derivatives, we can simply ignore their initialisation
1061 37 : const ADReal p = pressure;
1062 37 : const ADReal T = temperature;
1063 37 : const ADReal X = Xnacl;
1064 :
1065 : // FluidStateProperties data structure
1066 37 : std::vector<FluidStateProperties> fsp(_num_phases, FluidStateProperties(_num_components));
1067 37 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
1068 37 : FluidStateProperties & gas = fsp[_gas_phase_number];
1069 :
1070 : // Calculate equilibrium mass fractions in the two-phase state
1071 : ADReal Xco2, Yh2o;
1072 37 : equilibriumMassFractions(p, T, X, Xco2, Yh2o);
1073 :
1074 : // Save the mass fractions in the FluidStateMassFractions object
1075 37 : const ADReal Yco2 = 1.0 - Yh2o;
1076 74 : liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - Xco2;
1077 37 : liquid.mass_fraction[_gas_fluid_component] = Xco2;
1078 37 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
1079 37 : gas.mass_fraction[_gas_fluid_component] = Yco2;
1080 :
1081 : // Gas properties
1082 37 : gasProperties(pressure, temperature, fsp);
1083 :
1084 : // Liquid properties
1085 37 : const ADReal liquid_saturation = 1.0 - saturation;
1086 37 : const ADReal liquid_pressure = p - _pc.capillaryPressure(liquid_saturation, qp);
1087 37 : liquidProperties(liquid_pressure, T, X, fsp);
1088 :
1089 : // The total mass fraction of ncg (z) can now be calculated
1090 37 : const ADReal Z = (saturation * gas.density * Yco2 + liquid_saturation * liquid.density * Xco2) /
1091 74 : (saturation * gas.density + liquid_saturation * liquid.density);
1092 :
1093 37 : return Z.value();
1094 37 : }
1095 :
1096 : ADReal
1097 24 : PorousFlowBrineCO2::henryConstant(const ADReal & temperature, const ADReal & Xnacl) const
1098 : {
1099 : using std::pow;
1100 :
1101 : // Henry's constant for dissolution in water
1102 24 : const ADReal Kh_h2o = _brine_fp.henryConstant(temperature, _co2_henry);
1103 :
1104 : // The correction to salt is obtained through the salting out coefficient
1105 24 : const std::vector<Real> b{1.19784e-1, -7.17823e-4, 4.93854e-6, -1.03826e-8, 1.08233e-11};
1106 :
1107 : // Need temperature in Celsius
1108 : const ADReal Tc = temperature - _T_c2k;
1109 :
1110 24 : ADReal kb = 0.0;
1111 144 : for (unsigned int i = 0; i < b.size(); ++i)
1112 240 : kb += b[i] * pow(Tc, i);
1113 :
1114 : // Need salt mass fraction in molality
1115 48 : const ADReal xmol = Xnacl / (1.0 - Xnacl) / _Mnacl;
1116 :
1117 : // Henry's constant and its derivative wrt temperature and salt mass fraction
1118 48 : return Kh_h2o * pow(10.0, xmol * kb);
1119 24 : }
1120 :
1121 : ADReal
1122 6 : PorousFlowBrineCO2::enthalpyOfDissolutionGas(const ADReal & temperature, const ADReal & Xnacl) const
1123 : {
1124 : // Henry's constant
1125 6 : const ADReal Kh = henryConstant(temperature, Xnacl);
1126 :
1127 12 : ADReal hdis = -_R * temperature * temperature * Kh.derivatives()[_Tidx] / Kh / _Mco2;
1128 :
1129 : // Derivative of enthalpy of dissolution wrt temperature and xnacl requires the second
1130 : // derivatives of Henry's constant. For simplicity, approximate these numerically
1131 6 : const Real dT = temperature.value() * 1.0e-8;
1132 : const ADReal T2 = temperature + dT;
1133 6 : const ADReal Kh2 = henryConstant(T2, Xnacl);
1134 :
1135 : const Real dhdis_dT =
1136 18 : (-_R * T2 * T2 * Kh2.derivatives()[_Tidx] / Kh2 / _Mco2 - hdis).value() / dT;
1137 :
1138 6 : const Real dX = Xnacl.value() * 1.0e-8;
1139 : const ADReal X2 = Xnacl + dX;
1140 6 : const ADReal Kh3 = henryConstant(temperature, X2);
1141 :
1142 : const Real dhdis_dX =
1143 12 : (-_R * temperature * temperature * Kh3.derivatives()[_Tidx] / Kh3 / _Mco2 - hdis).value() /
1144 6 : dX;
1145 :
1146 12 : hdis.derivatives() = temperature.derivatives() * dhdis_dT + Xnacl.derivatives() * dhdis_dX;
1147 :
1148 6 : return hdis;
1149 : }
1150 :
1151 : ADReal
1152 627924 : PorousFlowBrineCO2::enthalpyOfDissolution(const ADReal & temperature) const
1153 : {
1154 : // Linear fit to model of Duan and Sun (2003) (in kJ/mol)
1155 1883772 : const ADReal delta_h = -58.3533 + 0.134519 * temperature;
1156 :
1157 : // Convert to J/kg
1158 627924 : return delta_h * 1000.0 / _Mco2;
1159 : }
1160 :
1161 : void
1162 4 : PorousFlowBrineCO2::smoothCubicInterpolation(
1163 : Real temperature, Real f0, Real df0, Real f1, Real df1, Real & value, Real & deriv) const
1164 : {
1165 : // Coefficients of cubic polynomial
1166 4 : const Real dT = _Tupper - _Tlower;
1167 :
1168 : const Real a = f0;
1169 4 : const Real b = df0 * dT;
1170 4 : const Real c = 3.0 * (f1 - f0) - (2.0 * df0 + df1) * dT;
1171 4 : const Real d = 2.0 * (f0 - f1) + (df0 + df1) * dT;
1172 :
1173 4 : const Real t2 = temperature * temperature;
1174 4 : const Real t3 = temperature * t2;
1175 :
1176 4 : value = a + b * temperature + c * t2 + d * t3;
1177 4 : deriv = b + 2.0 * c * temperature + 3.0 * d * t2;
1178 4 : }
1179 :
1180 : void
1181 636143 : PorousFlowBrineCO2::checkVariables(Real pressure, Real temperature) const
1182 : {
1183 : // The calculation of mass fractions is valid from 12C <= T <= 300C, and
1184 : // pressure less than 60 MPa
1185 636143 : if (temperature < 285.15 || temperature > 573.15)
1186 0 : mooseException(name() + ": temperature " + Moose::stringify(temperature) +
1187 : " is outside range 285.15 K <= T <= 573.15 K");
1188 :
1189 636143 : if (pressure > 6.0e7)
1190 0 : mooseException(name() + ": pressure " + Moose::stringify(pressure) +
1191 : " must be less than 60 MPa");
1192 636143 : }
|
