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 133865 : if (temperature.value() > 373.15)
404 0 : mooseError(name(),
405 : ": fugacityCoefficientsLowTemp() is not valid for T > 373.15K. Use "
406 : "fugacityCoefficientsHighTemp() instead");
407 :
408 : // Need pressure in bar
409 133865 : const ADReal pbar = pressure * 1.0e-5;
410 :
411 : // Molar volume in cm^3/mol
412 133865 : const ADReal V = _Mco2 / co2_density * 1.0e6;
413 :
414 : // Redlich-Kwong parameters
415 401595 : const ADReal aCO2 = 7.54e7 - 4.13e4 * temperature;
416 133865 : const Real bCO2 = 27.8;
417 133865 : const Real aCO2H2O = 7.89e7;
418 133865 : const Real bH2O = 18.18;
419 :
420 133865 : const ADReal t15 = std::pow(temperature, 1.5);
421 :
422 : // The fugacity coefficients for H2O and CO2
423 535460 : auto lnPhi = [V, aCO2, bCO2, t15, this](ADReal a, ADReal b)
424 : {
425 267730 : return std::log(V / (V - bCO2)) + b / (V - bCO2) -
426 535460 : 2.0 * a / (_Rbar * t15 * bCO2) * std::log((V + bCO2) / V) +
427 1338650 : aCO2 * b / (_Rbar * t15 * bCO2 * bCO2) * (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
428 133865 : };
429 :
430 401595 : const ADReal lnPhiH2O = lnPhi(aCO2H2O, bH2O) - std::log(pbar * V / (_Rbar * temperature));
431 133865 : const ADReal lnPhiCO2 = lnPhi(aCO2, bCO2) - std::log(pbar * V / (_Rbar * temperature));
432 :
433 133865 : fh2o = std::exp(lnPhiH2O);
434 133865 : fco2 = std::exp(lnPhiCO2);
435 133865 : }
436 :
437 : void
438 10402 : PorousFlowBrineCO2::fugacityCoefficientsHighTemp(const ADReal & pressure,
439 : const ADReal & temperature,
440 : const ADReal & co2_density,
441 : const ADReal & xco2,
442 : const ADReal & yh2o,
443 : ADReal & fco2,
444 : ADReal & fh2o) const
445 : {
446 10402 : if (temperature.value() <= 373.15)
447 0 : mooseError(name(),
448 : ": fugacityCoefficientsHighTemp() is not valid for T <= 373.15K. Use "
449 : "fugacityCoefficientsLowTemp() instead");
450 :
451 10402 : fh2o = fugacityCoefficientH2OHighTemp(pressure, temperature, co2_density, xco2, yh2o);
452 10402 : fco2 = fugacityCoefficientCO2HighTemp(pressure, temperature, co2_density, xco2, yh2o);
453 10402 : }
454 :
455 : ADReal
456 10403 : PorousFlowBrineCO2::fugacityCoefficientH2OHighTemp(const ADReal & pressure,
457 : const ADReal & temperature,
458 : const ADReal & co2_density,
459 : const ADReal & xco2,
460 : const ADReal & yh2o) const
461 : {
462 : // Need pressure in bar
463 10403 : const ADReal pbar = pressure * 1.0e-5;
464 : // Molar volume in cm^3/mol
465 10403 : const ADReal V = _Mco2 / co2_density * 1.0e6;
466 :
467 : // Redlich-Kwong parameters
468 10403 : const ADReal yco2 = 1.0 - yh2o;
469 10403 : const ADReal xh2o = 1.0 - xco2;
470 :
471 31209 : const ADReal aCO2 = 8.008e7 - 4.984e4 * temperature;
472 41612 : const ADReal aH2O = 1.337e8 - 1.4e4 * temperature;
473 10403 : const Real bCO2 = 28.25;
474 10403 : const Real bH2O = 15.7;
475 31209 : const ADReal KH2OCO2 = 1.427e-2 - 4.037e-4 * temperature;
476 31209 : const ADReal KCO2H2O = 0.4228 - 7.422e-4 * temperature;
477 10403 : const ADReal kH2OCO2 = KH2OCO2 * yh2o + KCO2H2O * yco2;
478 10403 : const ADReal kCO2H2O = KCO2H2O * yh2o + KH2OCO2 * yco2;
479 :
480 20806 : const ADReal aH2OCO2 = std::sqrt(aCO2 * aH2O) * (1.0 - kH2OCO2);
481 31209 : const ADReal aCO2H2O = std::sqrt(aCO2 * aH2O) * (1.0 - kCO2H2O);
482 :
483 10403 : const ADReal amix = yh2o * yh2o * aH2O + yh2o * yco2 * (aH2OCO2 + aCO2H2O) + yco2 * yco2 * aCO2;
484 10403 : const ADReal bmix = yh2o * bH2O + yco2 * bCO2;
485 :
486 10403 : const ADReal t15 = std::pow(temperature, 1.5);
487 :
488 0 : ADReal lnPhiH2O = bH2O / bmix * (pbar * V / (_Rbar * temperature) - 1.0) -
489 31209 : std::log(pbar * (V - bmix) / (_Rbar * temperature));
490 10403 : ADReal term3 = (2.0 * yh2o * aH2O + yco2 * (aH2OCO2 + aCO2H2O) -
491 10403 : yh2o * yco2 * std::sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O) * (yh2o - yco2) +
492 20806 : xh2o * xco2 * std::sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O)) /
493 : amix;
494 10403 : term3 -= bH2O / bmix;
495 20806 : term3 *= amix / (bmix * _Rbar * t15) * std::log(V / (V + bmix));
496 10403 : lnPhiH2O += term3;
497 :
498 10403 : return std::exp(lnPhiH2O);
499 : }
500 :
501 : ADReal
502 10403 : PorousFlowBrineCO2::fugacityCoefficientCO2HighTemp(const ADReal & pressure,
503 : const ADReal & temperature,
504 : const ADReal & co2_density,
505 : const ADReal & xco2,
506 : const ADReal & yh2o) const
507 : {
508 : // Need pressure in bar
509 10403 : const ADReal pbar = pressure * 1.0e-5;
510 : // Molar volume in cm^3/mol
511 10403 : const ADReal V = _Mco2 / co2_density * 1.0e6;
512 :
513 : // Redlich-Kwong parameters
514 10403 : const ADReal yco2 = 1.0 - yh2o;
515 10403 : const ADReal xh2o = 1.0 - xco2;
516 :
517 31209 : const ADReal aCO2 = 8.008e7 - 4.984e4 * temperature;
518 31209 : const ADReal aH2O = 1.337e8 - 1.4e4 * temperature;
519 10403 : const Real bCO2 = 28.25;
520 10403 : const Real bH2O = 15.7;
521 31209 : const ADReal KH2OCO2 = 1.427e-2 - 4.037e-4 * temperature;
522 31209 : const ADReal KCO2H2O = 0.4228 - 7.422e-4 * temperature;
523 10403 : const ADReal kH2OCO2 = KH2OCO2 * yh2o + KCO2H2O * yco2;
524 10403 : const ADReal kCO2H2O = KCO2H2O * yh2o + KH2OCO2 * yco2;
525 :
526 20806 : const ADReal aH2OCO2 = std::sqrt(aCO2 * aH2O) * (1.0 - kH2OCO2);
527 31209 : const ADReal aCO2H2O = std::sqrt(aCO2 * aH2O) * (1.0 - kCO2H2O);
528 :
529 10403 : const ADReal amix = yh2o * yh2o * aH2O + yh2o * yco2 * (aH2OCO2 + aCO2H2O) + yco2 * yco2 * aCO2;
530 10403 : const ADReal bmix = yh2o * bH2O + yco2 * bCO2;
531 :
532 10403 : const ADReal t15 = std::pow(temperature, 1.5);
533 :
534 0 : ADReal lnPhiCO2 = bCO2 / bmix * (pbar * V / (_Rbar * temperature) - 1.0) -
535 31209 : std::log(pbar * (V - bmix) / (_Rbar * temperature));
536 :
537 10403 : ADReal term3 = (2.0 * yco2 * aCO2 + yh2o * (aH2OCO2 + aCO2H2O) -
538 10403 : yh2o * yco2 * std::sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O) * (yh2o - yco2) +
539 20806 : xh2o * xco2 * std::sqrt(aH2O * aCO2) * (kCO2H2O - kH2OCO2)) /
540 : amix;
541 :
542 20806 : lnPhiCO2 += (term3 - bCO2 / bmix) * amix / (bmix * _Rbar * t15) * std::log(V / (V + bmix));
543 :
544 10403 : return std::exp(lnPhiCO2);
545 : }
546 :
547 : ADReal
548 10403 : PorousFlowBrineCO2::activityCoefficientH2O(const ADReal & temperature, const ADReal & xco2) const
549 : {
550 10403 : if (temperature.value() <= 373.15)
551 1 : return 1.0;
552 : else
553 : {
554 : const ADReal Tref = temperature - 373.15;
555 10402 : const ADReal xh2o = 1.0 - xco2;
556 20804 : const ADReal Am = -3.084e-2 * Tref + 1.927e-5 * Tref * Tref;
557 :
558 20804 : return std::exp((Am - 2.0 * Am * xh2o) * xco2 * xco2);
559 : }
560 : }
561 :
562 : ADReal
563 10403 : PorousFlowBrineCO2::activityCoefficientCO2(const ADReal & temperature, const ADReal & xco2) const
564 : {
565 10403 : if (temperature.value() <= 373.15)
566 1 : return 1.0;
567 : else
568 : {
569 : const ADReal Tref = temperature - 373.15;
570 10402 : const ADReal xh2o = 1.0 - xco2;
571 20804 : const ADReal Am = -3.084e-2 * Tref + 1.927e-5 * Tref * Tref;
572 :
573 20804 : return std::exp(2.0 * Am * xco2 * xh2o * xh2o);
574 : }
575 : }
576 :
577 : ADReal
578 133869 : PorousFlowBrineCO2::activityCoefficient(const ADReal & pressure,
579 : const ADReal & temperature,
580 : const ADReal & Xnacl) const
581 : {
582 : // Need pressure in bar
583 133869 : const ADReal pbar = pressure * 1.0e-5;
584 : // Need NaCl molality (mol/kg)
585 133869 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
586 :
587 401607 : const ADReal lambda = -0.411370585 + 6.07632013e-4 * temperature + 97.5347708 / temperature -
588 133869 : 0.0237622469 * pbar / temperature +
589 267738 : 0.0170656236 * pbar / (630.0 - temperature) +
590 133869 : 1.41335834e-5 * temperature * std::log(pbar);
591 :
592 267738 : const ADReal xi = 3.36389723e-4 - 1.9829898e-5 * temperature +
593 133869 : 2.12220830e-3 * pbar / temperature -
594 401607 : 5.24873303e-3 * pbar / (630.0 - temperature);
595 :
596 267738 : return std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
597 : }
598 :
599 : ADReal
600 10405 : PorousFlowBrineCO2::activityCoefficientHighTemp(const ADReal & temperature,
601 : const ADReal & Xnacl) const
602 : {
603 : // Need NaCl molality (mol/kg)
604 20810 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
605 :
606 : const ADReal T2 = temperature * temperature;
607 : const ADReal T3 = temperature * T2;
608 :
609 31215 : const ADReal lambda = 2.217e-4 * temperature + 1.074 / temperature + 2648.0 / T2;
610 41620 : const ADReal xi = 1.3e-5 * temperature - 20.12 / temperature + 5259.0 / T2;
611 :
612 41620 : return (1.0 - mnacl / _invMh2o) * std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
613 : }
614 :
615 : ADReal
616 144269 : PorousFlowBrineCO2::equilibriumConstantH2O(const ADReal & temperature) const
617 : {
618 : // Uses temperature in Celsius
619 : const ADReal Tc = temperature - _T_c2k;
620 : const ADReal Tc2 = Tc * Tc;
621 : const ADReal Tc3 = Tc2 * Tc;
622 : const ADReal Tc4 = Tc3 * Tc;
623 :
624 : ADReal logK0H2O;
625 :
626 144269 : if (Tc <= 99.0)
627 669330 : logK0H2O = -2.209 + 3.097e-2 * Tc - 1.098e-4 * Tc2 + 2.048e-7 * Tc3;
628 :
629 10403 : else if (Tc > 99.0 && Tc < 109.0)
630 : {
631 0 : const ADReal Tint = (Tc - 99.0) / 10.0;
632 : const ADReal Tint2 = Tint * Tint;
633 0 : logK0H2O = -0.0204026 + 0.0152513 * Tint + 0.417565 * Tint2 - 0.278636 * Tint * Tint2;
634 : }
635 :
636 : else // 109 <= Tc <= 300
637 62418 : logK0H2O = -2.1077 + 2.8127e-2 * Tc - 8.4298e-5 * Tc2 + 1.4969e-7 * Tc3 - 1.1812e-10 * Tc4;
638 :
639 144269 : return std::pow(10.0, logK0H2O);
640 : }
641 :
642 : ADReal
643 144269 : PorousFlowBrineCO2::equilibriumConstantCO2(const ADReal & temperature) const
644 : {
645 : // Uses temperature in Celsius
646 : const ADReal Tc = temperature - _T_c2k;
647 : const ADReal Tc2 = Tc * Tc;
648 : const ADReal Tc3 = Tc2 * Tc;
649 :
650 : ADReal logK0CO2;
651 :
652 144269 : if (Tc <= 99.0)
653 535464 : logK0CO2 = 1.189 + 1.304e-2 * Tc - 5.446e-5 * Tc2;
654 :
655 10403 : else if (Tc > 99.0 && Tc < 109.0)
656 : {
657 0 : const ADReal Tint = (Tc - 99.0) / 10.0;
658 : const ADReal Tint2 = Tint * Tint;
659 0 : logK0CO2 = 1.9462 + 2.25692e-2 * Tint - 9.49577e-3 * Tint2 - 6.77721e-3 * Tint * Tint2;
660 : }
661 :
662 : else // 109 <= Tc <= 300
663 52015 : logK0CO2 = 1.668 + 3.992e-3 * Tc - 1.156e-5 * Tc2 + 1.593e-9 * Tc3;
664 :
665 144269 : return std::pow(10.0, logK0CO2);
666 : }
667 :
668 : void
669 134276 : PorousFlowBrineCO2::equilibriumMoleFractions(const ADReal & pressure,
670 : const ADReal & temperature,
671 : const ADReal & Xnacl,
672 : ADReal & xco2,
673 : ADReal & yh2o) const
674 : {
675 134276 : if (temperature.value() <= _Tlower)
676 : {
677 133862 : equilibriumMoleFractionsLowTemp(pressure, temperature, Xnacl, xco2, yh2o);
678 : }
679 414 : else if (temperature.value() > _Tlower && temperature.value() < _Tupper)
680 : {
681 : // Cubic polynomial in this regime
682 2 : const Real Tint = (temperature.value() - _Tlower) / 10.0;
683 :
684 : // Equilibrium mole fractions and derivatives at the lower temperature
685 2 : ADReal Tlower = _Tlower;
686 2 : Moose::derivInsert(Tlower.derivatives(), _Tidx, 1.0);
687 :
688 : ADReal xco2_lower, yh2o_lower;
689 2 : equilibriumMoleFractionsLowTemp(pressure, Tlower, Xnacl, xco2_lower, yh2o_lower);
690 :
691 2 : const Real dxco2_dT_lower = xco2_lower.derivatives()[_Tidx];
692 2 : const Real dyh2o_dT_lower = yh2o_lower.derivatives()[_Tidx];
693 :
694 : // Equilibrium mole fractions and derivatives at the upper temperature
695 : Real xco2_upper, yh2o_upper;
696 2 : Real co2_density_upper = _co2_fp.rho_from_p_T(pressure.value(), _Tupper);
697 :
698 2 : solveEquilibriumMoleFractionHighTemp(
699 2 : pressure.value(), _Tupper, Xnacl.value(), co2_density_upper, xco2_upper, yh2o_upper);
700 :
701 : Real A, dA_dp, dA_dT, B, dB_dp, dB_dT, dB_dX;
702 2 : funcABHighTemp(pressure.value(),
703 2 : _Tupper,
704 : Xnacl.value(),
705 : co2_density_upper,
706 : xco2_upper,
707 : yh2o_upper,
708 : A,
709 : dA_dp,
710 : dA_dT,
711 : B,
712 : dB_dp,
713 : dB_dT,
714 : dB_dX);
715 :
716 2 : const Real dyh2o_dT_upper =
717 2 : ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
718 2 : const Real dxco2_dT_upper = dB_dT * (1.0 - yh2o_upper) - B * dyh2o_dT_upper;
719 :
720 : // The mole fractions in this regime are then found by interpolation
721 : Real xco2r, yh2or, dxco2_dT, dyh2o_dT;
722 2 : smoothCubicInterpolation(
723 : Tint, xco2_lower.value(), dxco2_dT_lower, xco2_upper, dxco2_dT_upper, xco2r, dxco2_dT);
724 2 : smoothCubicInterpolation(
725 : Tint, yh2o_lower.value(), dyh2o_dT_lower, yh2o_upper, dyh2o_dT_upper, yh2or, dyh2o_dT);
726 :
727 2 : xco2 = xco2r;
728 4 : Moose::derivInsert(xco2.derivatives(), _pidx, xco2_lower.derivatives()[_pidx]);
729 2 : Moose::derivInsert(xco2.derivatives(), _Tidx, dxco2_dT);
730 4 : Moose::derivInsert(xco2.derivatives(), _Xidx, xco2_lower.derivatives()[_Xidx]);
731 :
732 2 : yh2o = yh2or;
733 4 : Moose::derivInsert(yh2o.derivatives(), _pidx, yh2o_lower.derivatives()[_pidx]);
734 2 : Moose::derivInsert(yh2o.derivatives(), _Tidx, dyh2o_dT);
735 4 : Moose::derivInsert(yh2o.derivatives(), _Xidx, yh2o_lower.derivatives()[_Xidx]);
736 : }
737 : else
738 : {
739 : // CO2 density and derivatives wrt pressure and temperature
740 412 : const Real co2_density = _co2_fp.rho_from_p_T(pressure.value(), temperature.value());
741 :
742 : // Equilibrium mole fractions solved using iteration in this regime
743 : Real xco2r, yh2or;
744 412 : solveEquilibriumMoleFractionHighTemp(
745 : pressure.value(), temperature.value(), Xnacl.value(), co2_density, xco2r, yh2or);
746 :
747 : // Can use these in funcABHighTemp() to get derivatives analytically rather than by iteration
748 : Real A, dA_dp, dA_dT, B, dB_dp, dB_dT, dB_dX;
749 412 : funcABHighTemp(pressure.value(),
750 : temperature.value(),
751 : Xnacl.value(),
752 : co2_density,
753 : xco2r,
754 : yh2or,
755 : A,
756 : dA_dp,
757 : dA_dT,
758 : B,
759 : dB_dp,
760 : dB_dT,
761 : dB_dX);
762 :
763 412 : const Real dyh2o_dp =
764 412 : ((1.0 - B) * dA_dp + (A - 1.0) * A * dB_dp) / (1.0 - A * B) / (1.0 - A * B);
765 412 : const Real dxco2_dp = dB_dp * (1.0 - yh2or) - B * dyh2o_dp;
766 :
767 412 : const Real dyh2o_dT =
768 412 : ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
769 412 : const Real dxco2_dT = dB_dT * (1.0 - yh2or) - B * dyh2o_dT;
770 :
771 412 : const Real dyh2o_dX = ((A - 1.0) * A * dB_dX) / (1.0 - A * B) / (1.0 - A * B);
772 412 : const Real dxco2_dX = dB_dX * (1.0 - yh2or) - B * dyh2o_dX;
773 :
774 412 : xco2 = xco2r;
775 412 : Moose::derivInsert(xco2.derivatives(), _pidx, dxco2_dp);
776 412 : Moose::derivInsert(xco2.derivatives(), _Tidx, dxco2_dT);
777 412 : Moose::derivInsert(xco2.derivatives(), _Xidx, dxco2_dX);
778 :
779 412 : yh2o = yh2or;
780 412 : Moose::derivInsert(yh2o.derivatives(), _pidx, dyh2o_dp);
781 412 : Moose::derivInsert(yh2o.derivatives(), _Tidx, dyh2o_dT);
782 412 : Moose::derivInsert(yh2o.derivatives(), _Xidx, dyh2o_dX);
783 : }
784 134276 : }
785 :
786 : void
787 133864 : PorousFlowBrineCO2::equilibriumMoleFractionsLowTemp(const ADReal & pressure,
788 : const ADReal & temperature,
789 : const ADReal & Xnacl,
790 : ADReal & xco2,
791 : ADReal & yh2o) const
792 : {
793 133864 : if (temperature.value() > 373.15)
794 0 : mooseError(name(),
795 : ": equilibriumMoleFractionsLowTemp() is not valid for T > 373.15K. Use "
796 : "equilibriumMoleFractions() instead");
797 :
798 : // CO2 density and derivatives wrt pressure and temperature
799 133864 : const ADReal co2_density = _co2_fp.rho_from_p_T(pressure, temperature);
800 :
801 : // Assume infinite dilution (yh20 = 0 and xco2 = 0) in low temperature regime
802 : ADReal A, B;
803 133864 : funcABLowTemp(pressure, temperature, co2_density, A, B);
804 :
805 : // As the activity coefficient for CO2 in brine used in this regime isn't a 'true'
806 : // activity coefficient, we instead calculate the molality of CO2 in water, then
807 : // correct it for brine, and then calculate the mole fractions.
808 : // The mole fraction in pure water is
809 267728 : const ADReal yh2ow = (1.0 - B) / (1.0 / A - B);
810 133864 : const ADReal xco2w = B * (1.0 - yh2ow);
811 :
812 : // Molality of CO2 in pure water
813 267728 : const ADReal mco2w = xco2w * _invMh2o / (1.0 - xco2w);
814 : // Molality of CO2 in brine is then calculated using gamma
815 133864 : const ADReal gamma = activityCoefficient(pressure, temperature, Xnacl);
816 : const ADReal mco2 = mco2w / gamma;
817 :
818 : // Need NaCl molality (mol/kg)
819 133864 : const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
820 :
821 : // Mole fractions of CO2 and H2O in liquid and gas phases
822 267728 : const ADReal total_moles = 2.0 * mnacl + _invMh2o + mco2;
823 133864 : xco2 = mco2 / total_moles;
824 401592 : yh2o = A * (1.0 - xco2 - 2.0 * mnacl / total_moles);
825 133864 : }
826 :
827 : void
828 133864 : PorousFlowBrineCO2::funcABLowTemp(const ADReal & pressure,
829 : const ADReal & temperature,
830 : const ADReal & co2_density,
831 : ADReal & A,
832 : ADReal & B) const
833 : {
834 133864 : if (temperature.value() > 373.15)
835 0 : mooseError(name(),
836 : ": funcABLowTemp() is not valid for T > 373.15K. Use funcABHighTemp() instead");
837 :
838 : // Pressure in bar
839 133864 : const ADReal pbar = pressure * 1.0e-5;
840 :
841 : // Reference pressure and partial molar volumes
842 : const Real pref = 1.0;
843 133864 : const Real vCO2 = 32.6;
844 133864 : const Real vH2O = 18.1;
845 :
846 : const ADReal delta_pbar = pbar - pref;
847 133864 : const ADReal Rt = _Rbar * temperature;
848 :
849 : // Equilibrium constants
850 133864 : const ADReal K0H2O = equilibriumConstantH2O(temperature);
851 133864 : const ADReal K0CO2 = equilibriumConstantCO2(temperature);
852 :
853 : // Fugacity coefficients
854 : ADReal phiH2O, phiCO2;
855 133864 : fugacityCoefficientsLowTemp(pressure, temperature, co2_density, phiCO2, phiH2O);
856 :
857 267728 : A = K0H2O / (phiH2O * pbar) * std::exp(delta_pbar * vH2O / Rt);
858 401592 : B = phiCO2 * pbar / (_invMh2o * K0CO2) * std::exp(-delta_pbar * vCO2 / Rt);
859 133864 : }
860 :
861 : void
862 10401 : PorousFlowBrineCO2::funcABHighTemp(Real pressure,
863 : Real temperature,
864 : Real Xnacl,
865 : Real co2_density,
866 : Real xco2,
867 : Real yh2o,
868 : Real & A,
869 : Real & B) const
870 : {
871 10401 : if (temperature <= 373.15)
872 0 : mooseError(name(),
873 : ": funcABHighTemp() is not valid for T <= 373.15K. Use funcABLowTemp() instead");
874 :
875 : // Pressure in bar
876 10401 : const Real pbar = pressure * 1.0e-5;
877 : // Temperature in C
878 10401 : const Real Tc = temperature - _T_c2k;
879 :
880 : // Reference pressure and partial molar volumes
881 10401 : const Real pref = -1.9906e-1 + 2.0471e-3 * Tc + 1.0152e-4 * Tc * Tc - 1.4234e-6 * Tc * Tc * Tc +
882 10401 : 1.4168e-8 * Tc * Tc * Tc * Tc;
883 10401 : const Real vCO2 = 32.6 + 3.413e-2 * (Tc - 100.0);
884 10401 : const Real vH2O = 18.1 + 3.137e-2 * (Tc - 100.0);
885 :
886 10401 : const Real delta_pbar = pbar - pref;
887 10401 : const Real Rt = _Rbar * temperature;
888 :
889 : // Equilibrium constants
890 : // Use dummy ADReal temperature as derivatives aren't required
891 10401 : const ADReal T = temperature;
892 10401 : Real K0H2O = equilibriumConstantH2O(T).value();
893 10401 : Real K0CO2 = equilibriumConstantCO2(T).value();
894 :
895 : // Fugacity coefficients
896 : // Use dummy ADReal variables as derivatives aren't required
897 10401 : const ADReal p = pressure;
898 10401 : const ADReal rhoco2 = co2_density;
899 10401 : const ADReal x = xco2;
900 10401 : const ADReal y = yh2o;
901 :
902 : ADReal phiH2O, phiCO2;
903 10401 : fugacityCoefficientsHighTemp(p, T, rhoco2, x, y, phiCO2, phiH2O);
904 :
905 : // Activity coefficients
906 10401 : const Real gammaH2O = activityCoefficientH2O(T, x).value();
907 10401 : const Real gammaCO2 = activityCoefficientCO2(T, x).value();
908 :
909 : // Activity coefficient for CO2 in brine
910 : // Use dummy ADReal Xnacl as derivatives aren't required
911 10401 : const ADReal X = Xnacl;
912 10401 : const Real gamma = activityCoefficientHighTemp(T, X).value();
913 :
914 10401 : A = K0H2O * gammaH2O / (phiH2O.value() * pbar) * std::exp(delta_pbar * vH2O / Rt);
915 10401 : B = phiCO2.value() * pbar / (_invMh2o * K0CO2 * gamma * gammaCO2) *
916 10401 : std::exp(-delta_pbar * vCO2 / Rt);
917 10401 : }
918 :
919 : void
920 414 : PorousFlowBrineCO2::funcABHighTemp(Real pressure,
921 : Real temperature,
922 : Real Xnacl,
923 : Real co2_density,
924 : Real xco2,
925 : Real yh2o,
926 : Real & A,
927 : Real & dA_dp,
928 : Real & dA_dT,
929 : Real & B,
930 : Real & dB_dp,
931 : Real & dB_dT,
932 : Real & dB_dX) const
933 : {
934 414 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, xco2, yh2o, A, B);
935 :
936 : // Use finite differences for derivatives in the high temperature regime
937 : const Real dp = 1.0e-2;
938 : const Real dT = 1.0e-6;
939 : const Real dX = 1.0e-8;
940 :
941 : Real A2, B2;
942 414 : funcABHighTemp(pressure + dp, temperature, Xnacl, co2_density, xco2, yh2o, A2, B2);
943 414 : dA_dp = (A2 - A) / dp;
944 414 : dB_dp = (B2 - B) / dp;
945 :
946 414 : funcABHighTemp(pressure, temperature + dT, Xnacl, co2_density, xco2, yh2o, A2, B2);
947 414 : dA_dT = (A2 - A) / dT;
948 414 : dB_dT = (B2 - B) / dT;
949 :
950 414 : funcABHighTemp(pressure, temperature, Xnacl + dX, co2_density, xco2, yh2o, A2, B2);
951 414 : dB_dX = (B2 - B) / dX;
952 414 : }
953 :
954 : void
955 418 : PorousFlowBrineCO2::solveEquilibriumMoleFractionHighTemp(
956 : Real pressure, Real temperature, Real Xnacl, Real co2_density, Real & xco2, Real & yh2o) const
957 : {
958 : // Initial guess for yh2o and xco2 (from Spycher and Pruess (2010))
959 418 : Real y = _brine_fp.vaporPressure(temperature, 0.0) / pressure;
960 : Real x = 0.009;
961 :
962 : // Need salt mass fraction in molality
963 418 : const Real mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
964 :
965 : // If y > 1, then just use y = 1, x = 0 (only a gas phase)
966 418 : if (y >= 1.0)
967 : {
968 : y = 1.0;
969 : x = 0.0;
970 : }
971 : else
972 : {
973 : // Residual function for Newton-Raphson
974 8745 : auto fy = [mnacl, this](Real y, Real A, Real B) {
975 : return y -
976 8745 : (1.0 - B) * _invMh2o / ((1.0 / A - B) * (2.0 * mnacl + _invMh2o) + 2.0 * mnacl * B);
977 417 : };
978 :
979 : // Derivative of fy wrt y
980 4164 : auto dfy = [mnacl, this](Real A, Real B, Real dA, Real dB)
981 : {
982 4164 : const Real denominator = (1.0 / A - B) * (2.0 * mnacl + _invMh2o) + 2.0 * mnacl * B;
983 4164 : return 1.0 + _invMh2o * dB / denominator +
984 4164 : (1.0 - B) * _invMh2o *
985 4164 : (2.0 * mnacl * dB - (2.0 * mnacl + _invMh2o) * (dB + dA / A / A)) / denominator /
986 4164 : denominator;
987 417 : };
988 :
989 : Real A, B;
990 : Real dA, dB;
991 : const Real dy = 1.0e-8;
992 :
993 : // Solve for yh2o using Newton-Raphson method
994 : unsigned int iter = 0;
995 : const Real tol = 1.0e-12;
996 : const unsigned int max_its = 10;
997 417 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y, A, B);
998 :
999 4581 : while (std::abs(fy(y, A, B)) > tol)
1000 : {
1001 4164 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y, A, B);
1002 : // Finite difference derivatives of A and B wrt y
1003 4164 : funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y + dy, dA, dB);
1004 4164 : dA = (dA - A) / dy;
1005 4164 : dB = (dB - B) / dy;
1006 :
1007 4164 : y = y - fy(y, A, B) / dfy(A, B, dA, dB);
1008 :
1009 4164 : x = B * (1.0 - y);
1010 :
1011 : // Break if not converged and just use the value
1012 : if (iter > max_its)
1013 : break;
1014 : }
1015 : }
1016 :
1017 418 : yh2o = y;
1018 418 : xco2 = x;
1019 418 : }
1020 :
1021 : ADReal
1022 662712 : PorousFlowBrineCO2::partialDensityCO2(const ADReal & temperature) const
1023 : {
1024 : // This correlation uses temperature in C
1025 : const ADReal Tc = temperature - _T_c2k;
1026 : // The parial molar volume
1027 2650848 : const ADReal V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
1028 :
1029 1325424 : return 1.0e6 * _Mco2 / V;
1030 : }
1031 :
1032 : Real
1033 37 : PorousFlowBrineCO2::totalMassFraction(
1034 : Real pressure, Real temperature, Real Xnacl, Real saturation, unsigned int qp) const
1035 : {
1036 : // Check whether the input pressure and temperature are within the region of validity
1037 37 : checkVariables(pressure, temperature);
1038 :
1039 : // As we do not require derivatives, we can simply ignore their initialisation
1040 37 : const ADReal p = pressure;
1041 37 : const ADReal T = temperature;
1042 37 : const ADReal X = Xnacl;
1043 :
1044 : // FluidStateProperties data structure
1045 37 : std::vector<FluidStateProperties> fsp(_num_phases, FluidStateProperties(_num_components));
1046 37 : FluidStateProperties & liquid = fsp[_aqueous_phase_number];
1047 37 : FluidStateProperties & gas = fsp[_gas_phase_number];
1048 :
1049 : // Calculate equilibrium mass fractions in the two-phase state
1050 : ADReal Xco2, Yh2o;
1051 37 : equilibriumMassFractions(p, T, X, Xco2, Yh2o);
1052 :
1053 : // Save the mass fractions in the FluidStateMassFractions object
1054 37 : const ADReal Yco2 = 1.0 - Yh2o;
1055 74 : liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - Xco2;
1056 37 : liquid.mass_fraction[_gas_fluid_component] = Xco2;
1057 37 : gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
1058 37 : gas.mass_fraction[_gas_fluid_component] = Yco2;
1059 :
1060 : // Gas properties
1061 37 : gasProperties(pressure, temperature, fsp);
1062 :
1063 : // Liquid properties
1064 37 : const ADReal liquid_saturation = 1.0 - saturation;
1065 37 : const ADReal liquid_pressure = p - _pc.capillaryPressure(liquid_saturation, qp);
1066 37 : liquidProperties(liquid_pressure, T, X, fsp);
1067 :
1068 : // The total mass fraction of ncg (z) can now be calculated
1069 37 : const ADReal Z = (saturation * gas.density * Yco2 + liquid_saturation * liquid.density * Xco2) /
1070 74 : (saturation * gas.density + liquid_saturation * liquid.density);
1071 :
1072 37 : return Z.value();
1073 37 : }
1074 :
1075 : ADReal
1076 24 : PorousFlowBrineCO2::henryConstant(const ADReal & temperature, const ADReal & Xnacl) const
1077 : {
1078 : // Henry's constant for dissolution in water
1079 24 : const ADReal Kh_h2o = _brine_fp.henryConstant(temperature, _co2_henry);
1080 :
1081 : // The correction to salt is obtained through the salting out coefficient
1082 24 : const std::vector<Real> b{1.19784e-1, -7.17823e-4, 4.93854e-6, -1.03826e-8, 1.08233e-11};
1083 :
1084 : // Need temperature in Celsius
1085 : const ADReal Tc = temperature - _T_c2k;
1086 :
1087 24 : ADReal kb = 0.0;
1088 144 : for (unsigned int i = 0; i < b.size(); ++i)
1089 240 : kb += b[i] * std::pow(Tc, i);
1090 :
1091 : // Need salt mass fraction in molality
1092 48 : const ADReal xmol = Xnacl / (1.0 - Xnacl) / _Mnacl;
1093 :
1094 : // Henry's constant and its derivative wrt temperature and salt mass fraction
1095 48 : return Kh_h2o * std::pow(10.0, xmol * kb);
1096 24 : }
1097 :
1098 : ADReal
1099 6 : PorousFlowBrineCO2::enthalpyOfDissolutionGas(const ADReal & temperature, const ADReal & Xnacl) const
1100 : {
1101 : // Henry's constant
1102 6 : const ADReal Kh = henryConstant(temperature, Xnacl);
1103 :
1104 12 : ADReal hdis = -_R * temperature * temperature * Kh.derivatives()[_Tidx] / Kh / _Mco2;
1105 :
1106 : // Derivative of enthalpy of dissolution wrt temperature and xnacl requires the second
1107 : // derivatives of Henry's constant. For simplicity, approximate these numerically
1108 6 : const Real dT = temperature.value() * 1.0e-8;
1109 : const ADReal T2 = temperature + dT;
1110 6 : const ADReal Kh2 = henryConstant(T2, Xnacl);
1111 :
1112 : const Real dhdis_dT =
1113 18 : (-_R * T2 * T2 * Kh2.derivatives()[_Tidx] / Kh2 / _Mco2 - hdis).value() / dT;
1114 :
1115 6 : const Real dX = Xnacl.value() * 1.0e-8;
1116 : const ADReal X2 = Xnacl + dX;
1117 6 : const ADReal Kh3 = henryConstant(temperature, X2);
1118 :
1119 : const Real dhdis_dX =
1120 12 : (-_R * temperature * temperature * Kh3.derivatives()[_Tidx] / Kh3 / _Mco2 - hdis).value() /
1121 6 : dX;
1122 :
1123 12 : hdis.derivatives() = temperature.derivatives() * dhdis_dT + Xnacl.derivatives() * dhdis_dX;
1124 :
1125 6 : return hdis;
1126 : }
1127 :
1128 : ADReal
1129 627924 : PorousFlowBrineCO2::enthalpyOfDissolution(const ADReal & temperature) const
1130 : {
1131 : // Linear fit to model of Duan and Sun (2003) (in kJ/mol)
1132 1883772 : const ADReal delta_h = -58.3533 + 0.134519 * temperature;
1133 :
1134 : // Convert to J/kg
1135 627924 : return delta_h * 1000.0 / _Mco2;
1136 : }
1137 :
1138 : void
1139 4 : PorousFlowBrineCO2::smoothCubicInterpolation(
1140 : Real temperature, Real f0, Real df0, Real f1, Real df1, Real & value, Real & deriv) const
1141 : {
1142 : // Coefficients of cubic polynomial
1143 4 : const Real dT = _Tupper - _Tlower;
1144 :
1145 : const Real a = f0;
1146 4 : const Real b = df0 * dT;
1147 4 : const Real c = 3.0 * (f1 - f0) - (2.0 * df0 + df1) * dT;
1148 4 : const Real d = 2.0 * (f0 - f1) + (df0 + df1) * dT;
1149 :
1150 4 : const Real t2 = temperature * temperature;
1151 4 : const Real t3 = temperature * t2;
1152 :
1153 4 : value = a + b * temperature + c * t2 + d * t3;
1154 4 : deriv = b + 2.0 * c * temperature + 3.0 * d * t2;
1155 4 : }
1156 :
1157 : void
1158 636143 : PorousFlowBrineCO2::checkVariables(Real pressure, Real temperature) const
1159 : {
1160 : // The calculation of mass fractions is valid from 12C <= T <= 300C, and
1161 : // pressure less than 60 MPa
1162 636143 : if (temperature < 285.15 || temperature > 573.15)
1163 0 : mooseException(name() + ": temperature " + Moose::stringify(temperature) +
1164 : " is outside range 285.15 K <= T <= 573.15 K");
1165 :
1166 636143 : if (pressure > 6.0e7)
1167 0 : mooseException(name() + ": pressure " + Moose::stringify(pressure) +
1168 : " must be less than 60 MPa");
1169 636143 : }
|
