www.mooseframework.org
PorousFlowFluidStateBrineCO2.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 #include "BrineFluidProperties.h"
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<PorousFlowFluidStateFlashBase>();
18  params.addRequiredParam<UserObjectName>("brine_fp", "The name of the user object for brine");
19  params.addRequiredParam<UserObjectName>("co2_fp", "The name of the user object for CO2");
20  params.addCoupledVar("xnacl", 0, "The salt mass fraction in the brine (kg/kg)");
21  params.addClassDescription("Fluid state class for brine and CO2");
22  return params;
23 }
24 
26  : PorousFlowFluidStateFlashBase(parameters),
27 
28  _xnacl(_nodal_material ? coupledNodalValue("xnacl") : coupledValue("xnacl")),
29  _brine_fp(getUserObject<BrineFluidProperties>("brine_fp")),
30  _co2_fp(getUserObject<SinglePhaseFluidPropertiesPT>("co2_fp")),
31  _water_fp(_brine_fp.getComponent(BrineFluidProperties::WATER)),
32  _Mh2o(_brine_fp.molarMassH2O()),
33  _invMh2o(1.0 / _Mh2o),
34  _Mco2(_co2_fp.molarMass()),
35  _Mnacl(_brine_fp.molarMassNaCl()),
36  _Rbar(_R * 10.0)
37 {
38  // Check that the correct FluidProperties UserObjects have been provided
39  if (_brine_fp.fluidName() != "brine")
40  mooseError("Only a valid Brine FluidProperties UserObject can be provided in brine_fp");
41 
42  if (_co2_fp.fluidName() != "co2")
43  mooseError("Only a valid CO2 FluidProperties UserObject can be provided in co2_fp");
44 }
45 
46 void
48 {
49  // The FluidProperty objects use temperature in K
50  Real Tk = _temperature[_qp] + _T_c2k;
51  Real pressure = _gas_porepressure[_qp];
52  Real xnacl = _xnacl[_qp];
53 
54  // Mass fraction of CO2 in liquid and H2O in gas phases
55  Real X0, dX0_dp, dX0_dT, Y0, Y1, dY1_dp, dY1_dT;
56  massFractions(pressure, Tk, xnacl, X0, dX0_dp, dX0_dT, Y1, dY1_dp, dY1_dT);
57  Y0 = 1.0 - Y1;
58 
59  // Determine which phases are present based on the value of z
60  bool is_liquid = false;
61  bool is_gas = false;
62  bool is_twophase = false;
63 
64  if ((*_z[0])[_qp] <= X0)
65  {
66  // In this case, there is not enough CO2 to form a gas phase,
67  // so only a liquid phase is present
68  is_liquid = true;
69  }
70  else if ((*_z[0])[_qp] > X0 && (*_z[0])[_qp] < Y0)
71  {
72  // Two phases are present
73  is_twophase = true;
74  }
75  else // ((*_z[0])[_qp] >= Y0)
76  {
77  // In this case, there is not enough brine to form a liquid
78  // phase, so only a gas phase is present
79  is_gas = true;
80  }
81 
82  // Material must provide the following properties
83  Real vapor_mass_fraction = 0.0;
84  Real gas_saturation, liquid_saturation;
85  Real gas_density, liquid_density;
86  Real gas_viscosity, liquid_viscosity;
87 
88  // And the following derivatives
89  Real dX0_dz = 0.0, dY0_dp = 0.0, dY0_dT = 0.0, dY0_dz = 0.0;
90  Real dgas_density_dp = 0.0, dgas_density_dT = 0.0, dgas_density_dz = 0.0;
91  Real dliquid_density_dp = 0.0, dliquid_density_dT = 0.0, dliquid_density_dz = 0.0;
92  Real dgas_viscosity_dp = 0.0, dgas_viscosity_dT = 0.0, dgas_viscosity_dz = 0.0;
93  Real dliquid_viscosity_dp = 0.0, dliquid_viscosity_dT = 0.0, dliquid_viscosity_dz = 0.0;
94 
95  if (is_liquid)
96  {
97  X0 = (*_z[0])[_qp];
98  Y0 = 0.0;
99  gas_saturation = 0.0;
100  dX0_dp = 0.0;
101  dX0_dT = 0.0;
102  dX0_dz = 1.0;
103  gas_density = 0.0;
104  gas_viscosity = 1.0; // To guard against division by 0
105  }
106  else if (is_twophase)
107  {
108  // Two phase are present. Set mass equilibrium constants used in the
109  // calculation of vapor mass fraction
110  Y0 = 1.0 - Y1;
111  std::vector<Real> Ki(_num_components);
112  Ki[0] = Y0 / X0;
113  Ki[1] = (1.0 - Y0) / (1.0 - X0);
114  vapor_mass_fraction = vaporMassFraction(Ki);
115 
116  // Derivatives of mass fractions wrt PorousFlow variables
117  dY0_dp = -dY1_dp;
118  dY0_dT = -dY1_dT;
119  }
120  else // if (is_gas)
121  {
122  X0 = 0.0;
123  dX0_dp = 0.0;
124  dX0_dT = 0.0;
125  Y0 = (*_z[0])[_qp];
126  dY0_dz = 1.0;
127  vapor_mass_fraction = 1.0;
128  gas_saturation = 1.0;
129  liquid_density = 0.0;
130  liquid_viscosity = 1.0; // To guard against division by 0
131  }
132 
133  // Update all remaining mass fractions and derivatives
134  // Mass fraction of CO2 in gas and H2O in liquid phases
135  Real X1 = 1.0 - X0;
136  Y1 = 1.0 - Y0;
137 
138  // Calculate the gas density and viscosity in the single phase gas or two
139  // phase region
140  if (is_gas || is_twophase)
141  {
142  Real co2_density, dco2_density_dp, dco2_density_dT;
143  _co2_fp.rho_dpT(pressure, Tk, co2_density, dco2_density_dp, dco2_density_dT);
144 
145  // Gas density is given by the CO2 density - no correction due to the small amount of
146  // brine vapor is made
147  gas_density = co2_density;
148  dgas_density_dp = dco2_density_dp;
149  dgas_density_dT = dco2_density_dT;
150  dgas_density_dz = 0.0;
151 
152  Real co2_viscosity, dco2_viscosity_drho, dco2_viscosity_dT;
154  co2_density, Tk, dco2_density_dT, co2_viscosity, dco2_viscosity_drho, dco2_viscosity_dT);
155 
156  // Assume that the viscosity of the gas phase is a weighted sum of the
157  // individual viscosities
158  gas_viscosity = co2_viscosity;
159  dgas_viscosity_dp = dco2_viscosity_drho * dco2_density_dp;
160  dgas_viscosity_dT = dco2_viscosity_dT;
161  dgas_viscosity_dz = 0.0;
162  }
163 
164  // Calculate the saturation in the two phase case using the vapor mass
165  // fraction
166  if (is_twophase)
167  {
168  // Liquid density
169  Real co2_partial_density, dco2_partial_density_dT;
170  partialDensityCO2(Tk, co2_partial_density, dco2_partial_density_dT);
171  // Use old value of gas saturation to estimate liquid saturation
172  Real liqsat = 1.0;
173  if (!_is_initqp)
174  liqsat -= _saturation_old[_qp][_gas_phase_number];
175 
176  liquid_density =
177  1.0 / (X0 / co2_partial_density +
178  X1 / _brine_fp.rho(pressure - _pc_uo.capillaryPressure(liqsat), Tk, xnacl));
179 
180  // The gas saturation in the two phase case
181  gas_saturation = vapor_mass_fraction * liquid_density /
182  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
183  }
184 
185  // Calculate the saturations and pressures for each phase
186  liquid_saturation = 1.0 - gas_saturation;
187  Real liquid_porepressure = pressure - _pc_uo.capillaryPressure(liquid_saturation);
188 
189  // Calculate liquid density and viscosity if in the two phase or single phase
190  // liquid region, including a density correction due to the presence of dissolved
191  // CO2
192  if (is_liquid || is_twophase)
193  {
194  Real brine_density, dbrine_density_dp, dbrine_density_dT, dbrine_density_dx;
195  Real dliquid_viscosity_drho, dliquid_viscosity_dx;
196 
197  _brine_fp.rho_dpTx(liquid_porepressure,
198  Tk,
199  xnacl,
200  brine_density,
201  dbrine_density_dp,
202  dbrine_density_dT,
203  dbrine_density_dx);
204 
205  // The liquid density
206  Real co2_partial_density, dco2_partial_density_dT;
207  partialDensityCO2(Tk, co2_partial_density, dco2_partial_density_dT);
208  liquid_density = 1.0 / (X0 / co2_partial_density + X1 / brine_density);
209 
210  dliquid_density_dp =
211  (dX0_dp / brine_density + X1 * dbrine_density_dp / brine_density / brine_density -
212  dX0_dp / co2_partial_density) *
213  liquid_density * liquid_density;
214 
215  dliquid_density_dT =
216  (dX0_dT / brine_density + X1 * dbrine_density_dT / brine_density / brine_density -
217  dX0_dT / co2_partial_density +
218  X0 * dco2_partial_density_dT / co2_partial_density / co2_partial_density) *
219  liquid_density * liquid_density;
220 
221  dliquid_density_dz =
222  (dX0_dz / brine_density - dX0_dz / co2_partial_density) * liquid_density * liquid_density;
223 
224  // Liquid viscosity is just the brine viscosity
225  // Note: brine viscosity (and derivatives) requires water density (and derivatives)
226  Real water_density, dwater_density_dp, dwater_density_dT;
227  _water_fp.rho_dpT(liquid_porepressure, Tk, water_density, dwater_density_dp, dwater_density_dT);
228 
229  _brine_fp.mu_drhoTx(water_density,
230  Tk,
231  xnacl,
232  dwater_density_dT,
233  liquid_viscosity,
234  dliquid_viscosity_drho,
235  dliquid_viscosity_dT,
236  dliquid_viscosity_dx);
237 
238  // The derivative of viscosity wrt pressure is given by the chain rule
239  dliquid_viscosity_dp = dliquid_viscosity_drho * dwater_density_dp;
240  }
241 
242  // Save properties in FluidStateProperties vector
243  auto & liquid = _fsp[_aqueous_phase_number];
244  auto & gas = _fsp[_gas_phase_number];
245 
246  liquid.saturation = liquid_saturation;
247  gas.saturation = gas_saturation;
248  liquid.pressure = liquid_porepressure;
249  gas.pressure = _gas_porepressure[_qp];
250 
251  liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - X0;
252  liquid.mass_fraction[_gas_fluid_component] = X0;
253  gas.mass_fraction[_aqueous_fluid_component] = 1.0 - Y0;
254  gas.mass_fraction[_gas_fluid_component] = Y0;
255 
256  liquid.fluid_density = liquid_density;
257  gas.fluid_density = gas_density;
258  liquid.fluid_viscosity = liquid_viscosity;
259  gas.fluid_viscosity = gas_viscosity;
260 
261  // Derivatives wrt PorousFlow variables are required by the kernels
262  // Note: these don't need to be stateful so don't calculate them in
263  // initQpStatefulProperties
264  if (!_is_initqp)
265  {
266  // Derivative of gas saturation wrt variables
267  Real ds_dp = 0.0, ds_dT = 0.0, ds_dz = 0.0;
268  if (is_twophase)
269  {
270  Real K0 = Y0 / X0;
271  Real K1 = (1.0 - Y0) / (1.0 - X0);
272  Real dv_dz = (K1 - K0) / ((K0 - 1.0) * (K1 - 1.0));
273 
274  ds_dz = gas_density * liquid_density * dv_dz;
275  ds_dz /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
276  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
277 
278  Real dK0_dp = (-Y0 * dX0_dp - X0 * dY1_dp) / X0 / X0;
279  Real dK0_dT = (-Y0 * dX0_dT - X0 * dY1_dT) / X0 / X0;
280 
281  Real dK1_dp = (Y1 * dX0_dp + (1.0 - X0) * dY1_dp) / (1.0 - X0) / (1.0 - X0);
282  Real dK1_dT = (Y1 * dX0_dT + (1.0 - X0) * dY1_dT) / (1.0 - X0) / (1.0 - X0);
283 
284  Real dv_dp = (*_z[0])[_qp] * dK1_dp / (K1 - 1.0) / (K1 - 1.0) +
285  (1.0 - (*_z[0])[_qp]) * dK0_dp / (K0 - 1.0) / (K0 - 1.0);
286 
287  ds_dp = gas_density * liquid_density * dv_dp +
288  vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
289  (gas_density * dliquid_density_dp - dgas_density_dp * liquid_density);
290  ds_dp /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
291  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
292 
293  Real dv_dT = (*_z[0])[_qp] * dK1_dT / (K1 - 1.0) / (K1 - 1.0) +
294  (1.0 - (*_z[0])[_qp]) * dK0_dT / (K0 - 1.0) / (K0 - 1.0);
295 
296  ds_dT = gas_density * liquid_density * dv_dT +
297  vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
298  (gas_density * dliquid_density_dT - dgas_density_dT * liquid_density);
299  ds_dT /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
300  (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
301  }
302 
303  // Save derivatives in FluidStateProperties vector
304  liquid.dsaturation_dp = -ds_dp;
305  liquid.dsaturation_dT = -ds_dT;
306  liquid.dsaturation_dz = -ds_dz;
307  gas.dsaturation_dp = ds_dp;
308  gas.dsaturation_dT = ds_dT;
309  gas.dsaturation_dz = ds_dz;
310 
311  liquid.dmass_fraction_dp[_aqueous_fluid_component] = -dX0_dp;
312  liquid.dmass_fraction_dp[_gas_fluid_component] = dX0_dp;
313  liquid.dmass_fraction_dT[_aqueous_fluid_component] = -dX0_dT;
314  liquid.dmass_fraction_dT[_gas_fluid_component] = dX0_dT;
315  liquid.dmass_fraction_dz[_aqueous_fluid_component] = -dX0_dz;
316  liquid.dmass_fraction_dz[_gas_fluid_component] = dX0_dz;
317  gas.dmass_fraction_dp[_aqueous_fluid_component] = -dY0_dp;
318  gas.dmass_fraction_dp[_gas_fluid_component] = dY0_dp;
319  gas.dmass_fraction_dT[_aqueous_fluid_component] = -dY0_dT;
320  gas.dmass_fraction_dT[_gas_fluid_component] = dY0_dT;
321  gas.dmass_fraction_dz[_aqueous_fluid_component] = -dY0_dz;
322  gas.dmass_fraction_dz[_gas_fluid_component] = dY0_dz;
323 
324  liquid.dfluid_density_dp = dliquid_density_dp;
325  liquid.dfluid_density_dT = dliquid_density_dT;
326  liquid.dfluid_density_dz = dliquid_density_dz;
327  gas.dfluid_density_dp = dgas_density_dp;
328  gas.dfluid_density_dT = dgas_density_dT;
329  gas.dfluid_density_dz = dgas_density_dz;
330 
331  liquid.dfluid_viscosity_dp = dliquid_viscosity_dp;
332  liquid.dfluid_viscosity_dT = dliquid_viscosity_dT;
333  liquid.dfluid_viscosity_dz = dliquid_viscosity_dz;
334  gas.dfluid_viscosity_dp = dgas_viscosity_dp;
335  gas.dfluid_viscosity_dT = dgas_viscosity_dT;
336  gas.dfluid_viscosity_dz = dgas_viscosity_dz;
337  }
338 }
339 
340 void
342  Real temperature,
343  Real xnacl,
344  Real & xco2l,
345  Real & dxco2l_dp,
346  Real & dxco2l_dT,
347  Real & xh2og,
348  Real & dxh2og_dp,
349  Real & dxh2og_dT) const
350 {
351  // Pressure in bar
352  Real pbar = pressure * 1.0e-5;
353  // Pressure minus 1 bar
354  Real delta_pbar = pbar - 1.0;
355 
356  // Average partial molar volumes (cm^3/mol) as given by Sypcher, Pruess and Ennis-King (2003)
357  const Real vCO2 = 32.6;
358  const Real vH2O = 18.1;
359 
360  // NaCl molality (mol/kg)
361  Real mnacl = xnacl / (1.0 - xnacl) / _Mnacl;
362 
363  // Equilibrium constants
364  Real K0H2O, dK0H2O_dT, K0CO2, dK0CO2_dT;
365  equilibriumConstantH2O(temperature, K0H2O, dK0H2O_dT);
366  equilibriumConstantCO2(temperature, K0CO2, dK0CO2_dT);
367 
368  // Fugacity coefficients
369  Real phiH2O, dphiH2O_dp, dphiH2O_dT;
370  Real phiCO2, dphiCO2_dp, dphiCO2_dT;
371  fugacityCoefficientCO2(pressure, temperature, phiCO2, dphiCO2_dp, dphiCO2_dT);
372  fugacityCoefficientH2O(pressure, temperature, phiH2O, dphiH2O_dp, dphiH2O_dT);
373 
374  // Activity coefficient
375  Real gamma, dgamma_dp, dgamma_dT;
376  activityCoefficient(pressure, temperature, xnacl, gamma, dgamma_dp, dgamma_dT);
377 
378  Real A = K0H2O / (phiH2O * pbar) * std::exp(delta_pbar * vH2O / (_Rbar * temperature));
379  Real B =
380  phiCO2 * pbar / (_invMh2o * K0CO2) * std::exp(-delta_pbar * vCO2 / (_Rbar * temperature));
381 
382  Real dA_dp = (-1.0e-5 * K0H2O / pbar + 1.0e-5 * vH2O * K0H2O / (_Rbar * temperature) -
383  K0H2O * dphiH2O_dp / phiH2O) *
384  std::exp(delta_pbar * vH2O / (_Rbar * temperature)) / (pbar * phiH2O);
385  Real dB_dp = (1.0e-5 * phiCO2 + pbar * dphiCO2_dp -
386  1.0e-5 * vCO2 * pbar * phiCO2 / (_Rbar * temperature)) *
387  std::exp(-delta_pbar * vCO2 / (_Rbar * temperature)) / (_invMh2o * K0CO2);
388 
389  Real dA_dT = (dK0H2O_dT - dphiH2O_dT * K0H2O / phiH2O -
390  delta_pbar * vH2O * K0H2O / (_Rbar * temperature * temperature)) *
391  std::exp(delta_pbar * vH2O / (_Rbar * temperature)) / (pbar * phiH2O);
392  Real dB_dT = (-pbar * phiCO2 * dK0CO2_dT / K0CO2 + pbar * dphiCO2_dT +
393  delta_pbar * vCO2 * pbar * phiCO2 / (_Rbar * temperature * temperature)) *
394  std::exp(-delta_pbar * vCO2 / (_Rbar * temperature)) / (_invMh2o * K0CO2);
395 
396  // The mole fraction of H2O in the CO2-rich gas phase is then
397  Real yH2O = (1.0 - B) / (1.0 / A - B);
398  // The mole fraction of CO2 in the H2O-rich liquid phase is then (note: no salinty effect)
399  Real xCO2 = B * (1.0 - yH2O);
400  // The molality of CO2 in the H2O-rich liquid phase is (note: no salinty effect)
401  Real mCO2 = xCO2 * _invMh2o / (1.0 - xCO2);
402 
403  // The molality of CO2 in brine is then given by
404  Real mCO2b = mCO2 / gamma;
405  // The mole fraction of CO2 in brine is then
406  Real xCO2b = mCO2b / (2.0 * mnacl + _invMh2o + mCO2b);
407  // The mole fraction of H2O in the CO2-rich gas phase corrected for NaCl mole fraction is
408  Real yH2Ob = A * (1.0 - xCO2b - 2.0 * mnacl / (2.0 * mnacl + _invMh2o + mCO2b));
409 
410  // Convert the mole fractions to mass fractions and then update referenced values
411  // The mass fraction of H2O in gas (assume no salt in gas phase)
412  xh2og = yH2Ob * _Mh2o / (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2);
413 
414  // The number of moles of CO2 in 1kg of H2O
415  Real nco2 = xCO2b * (2.0 * mnacl + _invMh2o) / (1.0 - xCO2b);
416 
417  // The mass fraction of CO2 in liquid
418  xco2l = nco2 * _Mco2 / (1.0 + mnacl * _Mnacl + nco2 * _Mco2);
419 
420  // The derivatives of the mass fractions wrt pressure
421  Real dyH2O_dp = ((1.0 - B) * dA_dp + (A - 1.0) * A * dB_dp) / (1.0 - A * B) / (1.0 - A * B);
422  Real dxCO2_dp = dB_dp * (1.0 - yH2O) - B * dyH2O_dp;
423 
424  Real dmCO2_dp = _invMh2o * dxCO2_dp / (1.0 - xCO2) / (1.0 - xCO2);
425  Real dmCO2b_dp = dmCO2_dp / gamma - mCO2 * dgamma_dp / gamma / gamma;
426  Real dxCO2b_dp = (2.0 * mnacl + _invMh2o) * dmCO2b_dp / (2.0 * mnacl + _invMh2o + mCO2b) /
427  (2.0 * mnacl + _invMh2o + mCO2b);
428 
429  Real dyH2Ob_dp = (1.0 - xCO2b - 2.0 * mnacl / (2.0 * mnacl + _invMh2o + mCO2b)) * dA_dp -
430  A * dxCO2b_dp +
431  2.0 * A * mnacl * dmCO2b_dp / (2.0 * mnacl + _invMh2o + mCO2b) /
432  (2.0 * mnacl + _invMh2o + mCO2b);
433 
434  dxh2og_dp = _Mco2 * _Mh2o * dyH2Ob_dp / (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2) /
435  (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2);
436 
437  Real dnco2_dp = dxCO2b_dp * (2.0 * mnacl + _invMh2o) / (1.0 - xCO2b) / (1.0 - xCO2b);
438 
439  dxco2l_dp = (1.0 + mnacl * _Mnacl) * _Mco2 * dnco2_dp / (1.0 + mnacl * _Mnacl + nco2 * _Mco2) /
440  (1.0 + mnacl * _Mnacl + nco2 * _Mco2);
441 
442  // The derivatives of the mass fractions wrt temperature
443  Real dyH2O_dT = ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
444  Real dxCO2_dT = dB_dT * (1.0 - yH2O) - B * dyH2O_dT;
445 
446  Real dmCO2_dT = _invMh2o * dxCO2_dT / (1.0 - xCO2) / (1.0 - xCO2);
447  Real dmCO2b_dT = dmCO2_dT / gamma - mCO2 * dgamma_dT / gamma / gamma;
448  Real dxCO2b_dT = (2.0 * mnacl + _invMh2o) * dmCO2b_dT / (2.0 * mnacl + _invMh2o + mCO2b) /
449  (2.0 * mnacl + _invMh2o + mCO2b);
450 
451  Real dyH2Ob_dT = (1.0 - xCO2b - 2.0 * mnacl / (2.0 * mnacl + _invMh2o + mCO2b)) * dA_dT -
452  A * dxCO2b_dT +
453  2.0 * A * mnacl * dmCO2b_dT / (2.0 * mnacl + _invMh2o + mCO2b) /
454  (2.0 * mnacl + _invMh2o + mCO2b);
455 
456  dxh2og_dT = _Mco2 * _Mh2o * dyH2Ob_dT / (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2) /
457  (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2);
458 
459  Real dnco2_dT = dxCO2b_dT * (2.0 * mnacl + _invMh2o) / (1.0 - xCO2b) / (1.0 - xCO2b);
460 
461  dxco2l_dT = (1.0 + mnacl * _Mnacl) * _Mco2 * dnco2_dT / (1.0 + mnacl * _Mnacl + nco2 * _Mco2) /
462  (1.0 + mnacl * _Mnacl + nco2 * _Mco2);
463 }
464 
465 void
467  Real pressure, Real temperature, Real & fco2, Real & dfco2_dp, Real & dfco2_dT) const
468 {
469  // Need pressure in bar
470  Real pbar = pressure * 1.0e-5;
471  // CO2 density and derivatives wrt pressure and temperature
472  Real gas_density, dgas_density_dp, dgas_density_dT;
473  _co2_fp.rho_dpT(pressure, temperature, gas_density, dgas_density_dp, dgas_density_dT);
474  // Molar volume in cm^3/mol
475  Real V = _Mco2 / gas_density * 1.0e6;
476 
477  // Redlich-Kwong parameters
478  Real aCO2 = 7.54e7 - 4.13e4 * temperature;
479  Real bCO2 = 27.8;
480 
481  Real term1 = std::log(V / (V - bCO2)) + bCO2 / (V - bCO2) -
482  2.0 * aCO2 / (_Rbar * std::pow(temperature, 1.5) * bCO2) * std::log((V + bCO2) / V) +
483  aCO2 / (_Rbar * std::pow(temperature, 1.5) * bCO2) *
484  (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
485 
486  Real lnPhiCO2 = term1 - std::log(pbar * V / (_Rbar * temperature));
487  fco2 = std::exp(lnPhiCO2);
488 
489  // The derivative of the fugacity coefficient wrt pressure
490  Real dV_dp = -_Mco2 / gas_density / gas_density * dgas_density_dp * 1.0e6;
491  Real dterm1_dV =
492  (bCO2 * (bCO2 - 2.0 * V) / (bCO2 - V) / (bCO2 - V) +
493  aCO2 * (bCO2 + 2.0 * V) / (_Rbar * std::pow(temperature, 1.5) * (bCO2 + V) * (bCO2 + V))) /
494  V;
495  dfco2_dp = (dterm1_dV * dV_dp - 1.0e-5 / pbar - 1.0 / V * dV_dp) * fco2;
496 
497  // The derivative of the fugacity coefficient wrt temperature
498  Real dV_dT = -_Mco2 / gas_density / gas_density * dgas_density_dT * 1.0e6;
499  Real dterm1_dT = 3.0 * aCO2 * _Rbar * std::sqrt(temperature) * bCO2 * std::log((V + bCO2) / V) /
500  (_Rbar * std::pow(temperature, 1.5) * bCO2) /
501  (_Rbar * std::pow(temperature, 1.5) * bCO2) +
502  8.26e4 / (_Rbar * std::pow(temperature, 1.5) * bCO2) * std::log((V + bCO2) / V) -
503  1.5 * aCO2 * _Rbar * std::sqrt(temperature) * bCO2 *
504  (std::log((V + bCO2) / V) - bCO2 / (V + bCO2)) /
505  (_Rbar * std::pow(temperature, 1.5) * bCO2) /
506  (_Rbar * std::pow(temperature, 1.5) * bCO2) -
507  4.13e4 / (_Rbar * std::pow(temperature, 1.5) * bCO2) *
508  (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
509  dfco2_dT = (dterm1_dT + dterm1_dV * dV_dT - dV_dT / V + 1.0 / temperature) * fco2;
510 }
511 
512 void
514  Real pressure, Real temperature, Real & fh2o, Real & dfh2o_dp, Real & dfh2o_dT) const
515 {
516  // Need pressure in bar
517  Real pbar = pressure * 1.0e-5;
518  // CO2 density and derivatives wrt pressure and temperature
519  Real gas_density, dgas_density_dp, dgas_density_dT;
520  _co2_fp.rho_dpT(pressure, temperature, gas_density, dgas_density_dp, dgas_density_dT);
521  // Molar volume in cm^3/mol
522  Real V = _Mco2 / gas_density * 1.0e6;
523 
524  // Redlich-Kwong parameters
525  Real aCO2 = 7.54e7 - 4.13e4 * temperature;
526  Real bCO2 = 27.8;
527  Real aCO2H2O = 7.89e7;
528  Real bH2O = 18.18;
529 
530  Real term1 =
531  std::log(V / (V - bCO2)) + bH2O / (V - bCO2) -
532  2.0 * aCO2H2O / (_Rbar * std::pow(temperature, 1.5) * bCO2) * std::log((V + bCO2) / V) +
533  aCO2 * bH2O / (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2) *
534  (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
535 
536  Real lnPhiH2O = term1 - std::log(pbar * V / (_Rbar * temperature));
537  fh2o = std::exp(lnPhiH2O);
538 
539  // The derivative of the fugacity coefficient wrt pressure
540  Real dV_dp = -_Mco2 / gas_density / gas_density * dgas_density_dp * 1.0e6;
541  Real dterm1_dV = ((bCO2 * bCO2 - (bCO2 + bH2O) * V) / (bCO2 - V) / (bCO2 - V) +
542  (2.0 * aCO2H2O * (bCO2 + V) - aCO2 * bH2O) /
543  (_Rbar * std::pow(temperature, 1.5) * (bCO2 + V) * (bCO2 + V))) /
544  V;
545  dfh2o_dp = (dterm1_dV * dV_dp - 1.0e-5 / pbar - dV_dp / V) * fh2o;
546 
547  // The derivative of the fugacity coefficient wrt temperature
548  Real dV_dT = -_Mco2 / gas_density / gas_density * dgas_density_dT * 1.0e6;
549  Real dterm1_dT = 3.0 * _Rbar * std::sqrt(temperature) * bCO2 * aCO2H2O *
550  std::log((V + bCO2) / V) / (_Rbar * std::pow(temperature, 1.5) * bCO2) /
551  (_Rbar * std::pow(temperature, 1.5) * bCO2) -
552  1.5 * aCO2 * bH2O * _Rbar * std::sqrt(temperature) * bCO2 * bCO2 *
553  (std::log((V + bCO2) / V) - bCO2 / (V + bCO2)) /
554  (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2) /
555  (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2) -
556  4.13e4 * bH2O * (std::log((V + bCO2) / V) - bCO2 / (V + bCO2)) /
557  (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2);
558  dfh2o_dT = (dterm1_dT + dterm1_dV * dV_dT - dV_dT / V + 1.0 / temperature) * fh2o;
559 }
560 
561 void
563  Real temperature,
564  Real xnacl,
565  Real & gamma,
566  Real & dgamma_dp,
567  Real & dgamma_dT) const
568 {
569  // Need pressure in bar
570  Real pbar = pressure * 1.0e-5;
571  // Need NaCl molality (mol/kg)
572  Real mnacl = xnacl / (1.0 - xnacl) / _Mnacl;
573 
574  Real lambda = -0.411370585 + 6.07632013e-4 * temperature + 97.5347708 / temperature -
575  0.0237622469 * pbar / temperature + 0.0170656236 * pbar / (630.0 - temperature) +
576  1.41335834e-5 * temperature * std::log(pbar);
577 
578  Real xi = 3.36389723e-4 - 1.9829898e-5 * temperature + 2.12220830e-3 * pbar / temperature -
579  5.24873303e-3 * pbar / (630.0 - temperature);
580 
581  gamma = std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
582 
583  // Derivative wrt pressure
584  Real dlambda_dp, dxi_dp;
585  dlambda_dp = -0.0237622469 / temperature + 0.0170656236 / (630.0 - temperature) +
586  1.41335834e-5 * temperature / pbar;
587  dxi_dp = 2.12220830e-3 / temperature - 5.24873303e-3 / (630.0 - temperature);
588  dgamma_dp = (2.0 * mnacl * dlambda_dp + mnacl * mnacl * dxi_dp) *
589  std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl) * 1.0e-5;
590 
591  // Derivative wrt temperature
592  Real dlambda_dT, dxi_dT;
593  dlambda_dT = 6.07632013e-4 - 97.5347708 / temperature / temperature +
594  0.0237622469 * pbar / temperature / temperature +
595  0.0170656236 * pbar / (630.0 - temperature) / (630.0 - temperature) +
596  1.41335834e-5 * std::log(pbar);
597  dxi_dT = -1.9829898e-5 - 2.12220830e-3 * pbar / temperature / temperature -
598  5.24873303e-3 * pbar / (630.0 - temperature) / (630.0 - temperature);
599  dgamma_dT = (2.0 * mnacl * dlambda_dT + mnacl * mnacl * dxi_dT) *
600  std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
601 }
602 
603 void
605  Real & kh2o,
606  Real & dkh2o_dT) const
607 {
608  // Uses temperature in Celcius
609  Real Tc = temperature - _T_c2k;
610 
611  Real logK0H2O = -2.209 + 3.097e-2 * Tc - 1.098e-4 * Tc * Tc + 2.048e-7 * Tc * Tc * Tc;
612  Real dlogK0H2O = 3.097e-2 - 2.196e-4 * Tc + 6.144e-7 * Tc * Tc;
613 
614  kh2o = std::pow(10.0, logK0H2O);
615  dkh2o_dT = std::log(10.0) * dlogK0H2O * kh2o;
616 }
617 
618 void
620  Real & kco2,
621  Real & dkco2_dT) const
622 {
623  // Uses temperature in Celcius
624  Real Tc = temperature - _T_c2k;
625 
626  Real logK0CO2 = 1.189 + 1.304e-2 * Tc - 5.446e-5 * Tc * Tc;
627  Real dlogK0CO2 = 1.304e-2 - 1.0892e-4 * Tc;
628 
629  kco2 = std::pow(10.0, logK0CO2);
630  dkco2_dT = std::log(10.0) * dlogK0CO2 * kco2;
631 }
632 
633 void
635  Real & partial_density,
636  Real & dpartial_density_dT) const
637 {
638  // This correlation uses temperature in C
639  Real Tc = temperature - _T_c2k;
640  // The parial molar volume
641  Real V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
642  Real dV_dT = -9.585e-2 + 1.748e-3 * Tc - 1.5132e-6 * Tc * Tc;
643 
644  partial_density = 1.0e6 * _Mco2 / V;
645  dpartial_density_dT = -1.0e6 * _Mco2 * dV_dT / V / V;
646 }
virtual std::string fluidName() const =0
Fluid name.
PorousFlowFluidStateBrineCO2(const InputParameters &parameters)
virtual void rho_dpTx(Real pressure, Real temperature, Real xnacl, Real &rho, Real &drho_dp, Real &drho_dT, Real &drho_dx) const override
Density and its derivatives wrt pressure, temperature and mass fraction.
const MaterialProperty< Real > & _temperature
Temperature.
void fugacityCoefficientCO2(Real pressure, Real temperature, Real &fco2, Real &dfco2_dp, Real &dfco2_dT) const
Fugacity coefficient for CO2.
virtual void rho_dpT(Real pressure, Real temperature, Real &rho, Real &drho_dp, Real &drho_dT) const =0
Density and its derivatives wrt pressure and temperature.
const Real _Mnacl
Molar mass of NaCL.
void fugacityCoefficientH2O(Real pressure, Real temperature, Real &fh2o, Real &dfh2o_dp, Real &dfh2o_dT) const
Fugacity coefficient for H2O.
const SinglePhaseFluidPropertiesPT & _water_fp
Fluid properties UserObject for H20.
const std::string temperature
Definition: NS.h:25
const BrineFluidProperties & _brine_fp
Fluid properties UserObject for brine.
virtual std::string fluidName() const override
Fluid name.
Brine (NaCl in H2O) fluid properties as a function of pressure (Pa), temperature (K) and NaCl mass fr...
const VariableValue & _gas_porepressure
Porepressure.
void activityCoefficient(Real pressure, Real temperature, Real xnacl, Real &gamma, Real &dgamma_dp, Real &dgamma_dT) const
Activity coefficient for CO2 in brine.
const MaterialProperty< std::vector< Real > > & _saturation_old
Old value of saturation.
void massFractions(Real pressure, Real temperature, Real xnacl, Real &xco2l, Real &dxco2l_dp, Real &dxco2l_dT, Real &xh2og, Real &dxh2og_dp, Real &dxh2og_dT) const
Mass fractions of CO2 and brine calculated using mutual solubilities given by Spycher, Pruess and Ennis-King, CO2-H2O mixtures in the geological sequestration of CO2.
Common class for single phase fluid properties using a pressure and temperature formulation.
virtual Real rho(Real pressure, Real temperature, Real xnacl) const override
Density.
const VariableValue & _xnacl
Salt mass fraction (kg/kg)
virtual Real capillaryPressure(Real saturation) const
Capillary pressure is calculated as a function of true saturation.
const unsigned int _gas_fluid_component
Fluid component number of the gas phase.
const PorousFlowCapillaryPressure & _pc_uo
Capillary pressure UserObject.
void partialDensityCO2(Real temperature, Real &partial_density, Real &dpartial_density_dT) const
Partial density of dissolved CO2 From Garcia, Density of aqueous solutions of CO2, LBNL-49023 (2001)
const Real _Mh2o
Molar mass of water (kg/mol)
std::vector< FluidStateProperties > _fsp
FluidStateProperties data structure.
const unsigned int _aqueous_phase_number
Phase number of the aqueous phase.
virtual void mu_drhoT_from_rho_T(Real density, Real temperature, Real ddensity_dT, Real &mu, Real &dmu_drho, Real &dmu_dT) const =0
Dynamic viscosity and its derivatives wrt density and temperature.
const SinglePhaseFluidPropertiesPT & _co2_fp
Fluid properties UserObject for CO2.
virtual void thermophysicalProperties() override
Calculates all required thermophysical properties and derivatives for each phase and fluid component...
const unsigned int _gas_phase_number
Phase number of the gas phase.
bool _is_initqp
Flag to indicate whether to calculate stateful properties.
InputParameters validParams< PorousFlowFluidStateBrineCO2 >()
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const Real _invMh2o
Inverse of molar mass of H2O (mol/kg)
const unsigned int _aqueous_fluid_component
Fluid component number of the aqueous component.
InputParameters validParams< PorousFlowFluidStateFlashBase >()
const Real _T_c2k
Conversion from degrees Celsius to degrees Kelvin.
const unsigned int _num_components
Number of components.
const std::string pressure
Definition: NS.h:24
const Real _Mco2
Molar mass of CO2 (kg/mol)
void equilibriumConstantH2O(Real temperature, Real &kh2o, Real &dkh2o_dT) const
Equilibrium constant for H2O from Spycher, Pruess and Ennis-King, CO2-H2O mixtures in the geological ...
const Real _Rbar
Molar gas constant in bar cm^3 /(K mol)
virtual void mu_drhoTx(Real water_density, Real temperature, Real xnacl, Real dwater_density_dT, Real &mu, Real &dmu_drho, Real &dmu_dT, Real &dmu_dx) const override
Dynamic viscosity and its derivatives wrt density, temperature and mass fraction. ...
void equilibriumConstantCO2(Real temperature, Real &kco2, Real &dkco2_dT) const
Equilibrium constant for CO2 from Spycher, Pruess and Ennis-King, CO2-H2O mixtures in the geological ...
std::vector< const VariableValue * > _z
Total mass fraction(s) of the gas component(s) summed over all phases.
Base class for fluid states using a persistent set of primary variables for the mutliphase, multicomponent case.
Real vaporMassFraction(std::vector< Real > &Ki) const
Solves Rachford-Rice equation to provide vapor mass fraction.