Specialized class for brine and CO2 including calculation of mutual solubility of the two fluids using a high-accuracy fugacity-based formulation. More...

#include <PorousFlowBrineCO2.h>

Inheritance diagram for PorousFlowBrineCO2:

[legend]

Public Types
typedef DataFileName	DataFileParameterType

Public Member Functions
	PorousFlowBrineCO2 (const InputParameters &parameters)

virtual std::string	fluidStateName () const override
	Name of FluidState. More...

void	thermophysicalProperties (Real pressure, Real temperature, Real Xnacl, Real Z, unsigned int qp, std::vector< FluidStateProperties > &fsp) const override
	Determines the complete thermophysical state of the system for a given set of primary variables. More...

void	thermophysicalProperties (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, const ADReal &Z, unsigned int qp, std::vector< FluidStateProperties > &fsp) const override

void	equilibriumMoleFractions (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, ADReal &xco2, ADReal &yh2o) const
	Mole fractions of CO2 in brine and water vapor in CO2 at equilibrium. More...

void	equilibriumMassFractions (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, ADReal &Xco2, ADReal &Yh2o) const
	Mass fractions of CO2 in brine and water vapor in CO2 at equilibrium. More...

void	massFractions (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, const ADReal &Z, FluidStatePhaseEnum &phase_state, std::vector< FluidStateProperties > &fsp) const
	Mass fractions of CO2 and H2O in both phases, as well as derivatives wrt PorousFlow variables. More...

void	gasProperties (const ADReal &pressure, const ADReal &temperature, std::vector< FluidStateProperties > &fsp) const
	Thermophysical properties of the gaseous state. More...

void	liquidProperties (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, std::vector< FluidStateProperties > &fsp) const
	Thermophysical properties of the liquid state. More...

ADReal	saturation (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, const ADReal &Z, std::vector< FluidStateProperties > &fsp) const
	Gas saturation in the two-phase region. More...

void	twoPhaseProperties (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, const ADReal &Z, unsigned int qp, std::vector< FluidStateProperties > &fsp) const
	Gas and liquid properties in the two-phase region. More...

void	fugacityCoefficientsLowTemp (const ADReal &pressure, const ADReal &temperature, const ADReal &co2_density, ADReal &fco2, ADReal &fh2o) const
	Fugacity coefficients for H2O and CO2 for T <= 100C Eq. More...

ADReal	activityCoefficientH2O (const ADReal &temperature, const ADReal &xco2) const
	Activity coefficient of H2O Eq. More...

ADReal	activityCoefficientCO2 (const ADReal &temperature, const ADReal &xco2) const
	Activity coefficient of CO2 Eq. More...

ADReal	activityCoefficient (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl) const
	Activity coefficient for CO2 in brine. More...

ADReal	activityCoefficientHighTemp (const ADReal &temperature, const ADReal &Xnacl) const
	Activity coefficient for CO2 in brine used in the elevated temperature formulation. More...

ADReal	equilibriumConstantH2O (const ADReal &temperature) const
	Equilibrium constant for H2O For temperatures 12C <= T <= 99C, uses Spycher, Pruess and Ennis-King (2003) For temperatures 109C <= T <= 300C, uses Spycher and Pruess (2010) For temperatures 99C < T < 109C, the value is calculated by smoothly interpolating the two formulations. More...

ADReal	equilibriumConstantCO2 (const ADReal &temperature) const
	Equilibrium constant for CO2 For temperatures 12C <= T <= 99C, uses Spycher, Pruess and Ennis-King (2003) For temperatures 109C <= T <= 300C, uses Spycher and Pruess (2010) For temperatures 99C < T < 109C, the value is calculated by smoothly interpolating the two formulations. More...

ADReal	partialDensityCO2 (const ADReal &temperature) const
	Partial density of dissolved CO2 From Garcia, Density of aqueous solutions of CO2, LBNL-49023 (2001) More...

virtual Real	totalMassFraction (Real pressure, Real temperature, Real Xnacl, Real saturation, unsigned int qp) const override
	Total mass fraction of fluid component summed over all phases in the two-phase state for a specified gas saturation. More...

unsigned int	saltComponentIndex () const
	The index of the salt component. More...

ADReal	henryConstant (const ADReal &temperature, const ADReal &Xnacl) const
	Henry's constant of dissolution of gas phase CO2 in brine. More...

ADReal	enthalpyOfDissolutionGas (const ADReal &temperature, const ADReal &Xnacl) const
	Enthalpy of dissolution of gas phase CO2 in brine calculated using Henry's constant From Himmelblau, Partial molal heats and entropies of solution for gases dissolved in water from the freezing to the near critical point, J. More...

ADReal	enthalpyOfDissolution (const ADReal &temperature) const
	Enthalpy of dissolution of CO2 in brine calculated using linear fit to model of Duan and Sun, An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar, Chemical geology, 193, 257–271 (2003). More...

void	equilibriumMoleFractionsLowTemp (const ADReal &pressure, const ADReal &temperature, const ADReal &Xnacl, ADReal &xco2, ADReal &yh2o) const
	Mole fractions of CO2 in brine and water vapor in CO2 at equilibrium in the low temperature regime (T <= 99C). More...

void	solveEquilibriumMoleFractionHighTemp (Real pressure, Real temperature, Real Xnacl, Real co2_density, Real &xco2, Real &yh2o) const
	Function to solve for yh2o and xco2 iteratively in the elevated temperature regime (T > 100C) More...

unsigned int	aqueousComponentIndex () const
	The index of the aqueous fluid component. More...

unsigned int	gasComponentIndex () const
	The index of the gas fluid component. More...

void	phaseState (Real Zi, Real Xi, Real Yi, FluidStatePhaseEnum &phase_state) const
	Determines the phase state gven the total mass fraction and equilibrium mass fractions. More...

unsigned int	getPressureIndex () const

unsigned int	getTemperatureIndex () const

unsigned int	getZIndex () const

unsigned int	getXIndex () const

Real	rachfordRice (Real vf, std::vector< Real > &Zi, std::vector< Real > &Ki) const
	Rachford-Rice equation for vapor fraction. More...

Real	rachfordRiceDeriv (Real vf, std::vector< Real > &Zi, std::vector< Real > &Ki) const
	Derivative of Rachford-Rice equation wrt vapor fraction. More...

Real	vaporMassFraction (Real Z0, Real K0, Real K1) const
	Solves Rachford-Rice equation to provide vapor mass fraction. More...

ADReal	vaporMassFraction (const ADReal &Z0, const ADReal &K0, const ADReal &K1) const

Real	vaporMassFraction (std::vector< Real > &Zi, std::vector< Real > &Ki) const

void	initialize () final

void	execute () final

void	finalize () final

unsigned int	numPhases () const
	The maximum number of phases in this model. More...

unsigned int	numComponents () const
	The maximum number of components in this model. More...

unsigned int	aqueousPhaseIndex () const
	The index of the aqueous phase. More...

unsigned int	gasPhaseIndex () const
	The index of the gas phase. More...

void	clearFluidStateProperties (std::vector< FluidStateProperties > &fsp) const
	Clears the contents of the FluidStateProperties data structure. More...

SubProblem &	getSubProblem () const

bool	shouldDuplicateInitialExecution () const

virtual Real	spatialValue (const Point &) const

virtual const std::vector< Point >	spatialPoints () const

void	gatherSum (T &value)

void	gatherMax (T &value)

void	gatherMin (T &value)

void	gatherProxyValueMax (T1 &proxy, T2 &value)

void	gatherProxyValueMin (T1 &proxy, T2 &value)

void	setPrimaryThreadCopy (UserObject *primary)

UserObject *	primaryThreadCopy ()

std::set< UserObjectName >	getDependObjects () const

virtual bool	needThreadedCopy () const

const std::set< std::string > &	getRequestedItems () override

const std::set< std::string > &	getSuppliedItems () override

unsigned int	systemNumber () const

virtual bool	enabled () const

std::shared_ptr< MooseObject >	getSharedPtr ()

std::shared_ptr< const MooseObject >	getSharedPtr () const

bool	isKokkosObject (IsKokkosObjectKey &&) const

MooseApp &	getMooseApp () const

const std::string &	type () const

const std::string &	name () const

std::string	typeAndName () const

MooseObjectParameterName	uniqueParameterName (const std::string &parameter_name) const

MooseObjectName	uniqueName () const

const InputParameters &	parameters () const

const hit::Node *	getHitNode () const

bool	hasBase () const

const std::string &	getBase () const

const T &	getParam (const std::string &name) const

std::vector< std::pair< T1, T2 > >	getParam (const std::string &param1, const std::string &param2) const

const T *	queryParam (const std::string &name) const

const T &	getRenamedParam (const std::string &old_name, const std::string &new_name) const

T	getCheckedPointerParam (const std::string &name, const std::string &error_string="") const

bool	isParamValid (const std::string &name) const

bool	isParamSetByUser (const std::string &name) const

void	connectControllableParams (const std::string &parameter, const std::string &object_type, const std::string &object_name, const std::string &object_parameter) const

void	paramError (const std::string &param, Args... args) const

void	paramWarning (const std::string &param, Args... args) const

void	paramWarning (const std::string &param, Args... args) const

void	paramInfo (const std::string &param, Args... args) const

std::string	messagePrefix (const bool hit_prefix=true) const

std::string	errorPrefix (const std::string &) const

void	mooseError (Args &&... args) const

void	mooseDocumentedError (const std::string &repo_name, const unsigned int issue_num, Args &&... args) const

void	mooseErrorNonPrefixed (Args &&... args) const

void	mooseWarning (Args &&... args) const

void	mooseWarning (Args &&... args) const

void	mooseWarningNonPrefixed (Args &&... args) const

void	mooseWarningNonPrefixed (Args &&... args) const

void	mooseDeprecated (Args &&... args) const

void	mooseDeprecated (Args &&... args) const

void	mooseInfo (Args &&... args) const

void	callMooseError (std::string msg, const bool with_prefix, const hit::Node *node=nullptr) const

std::string	getDataFileName (const std::string &param) const

std::string	getDataFileNameByName (const std::string &relative_path) const

std::string	getDataFilePath (const std::string &relative_path) const

virtual void	initialSetup ()

virtual void	timestepSetup ()

virtual void	jacobianSetup ()

virtual void	residualSetup ()

virtual void	customSetup (const ExecFlagType &)

const ExecFlagEnum &	getExecuteOnEnum () const

UserObjectName	getUserObjectName (const std::string &param_name) const

const T &	getUserObject (const std::string &param_name, bool is_dependency=true) const

const T &	getUserObjectByName (const UserObjectName &object_name, bool is_dependency=true) const

const UserObject &	getUserObjectBase (const std::string &param_name, bool is_dependency=true) const

const UserObject &	getUserObjectBaseByName (const UserObjectName &object_name, bool is_dependency=true) const

const std::vector< MooseVariableScalar *> &	getCoupledMooseScalarVars ()

const std::set< TagID > &	getScalarVariableCoupleableVectorTags () const

const std::set< TagID > &	getScalarVariableCoupleableMatrixTags () const

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialProperty (const std::string &name, MaterialData &material_data, const unsigned int state=0)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialProperty (const std::string &name, const unsigned int state=0)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialProperty (const std::string &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialProperty (const std::string &name, MaterialData &material_data, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialProperty (const std::string &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialProperty (const std::string &name, const unsigned int state=0)

const ADMaterialProperty< T > &	getADMaterialProperty (const std::string &name, MaterialData &material_data)

const ADMaterialProperty< T > &	getADMaterialProperty (const std::string &name)

const ADMaterialProperty< T > &	getADMaterialProperty (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOld (const std::string &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOld (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOld (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOlder (const std::string &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOlder (const std::string &name)

const MaterialProperty< T > &	getMaterialPropertyOlder (const std::string &name)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialPropertyByName (const MaterialPropertyName &name, MaterialData &material_data, const unsigned int state)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const GenericMaterialProperty< T, is_ad > &	getGenericMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialPropertyByName (const MaterialPropertyName &name, MaterialData &material_data, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const MaterialProperty< T > &	getMaterialPropertyByName (const MaterialPropertyName &name, const unsigned int state=0)

const ADMaterialProperty< T > &	getADMaterialPropertyByName (const MaterialPropertyName &name, MaterialData &material_data)

const ADMaterialProperty< T > &	getADMaterialPropertyByName (const MaterialPropertyName &name)

const ADMaterialProperty< T > &	getADMaterialPropertyByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOldByName (const MaterialPropertyName &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOldByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOldByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOlderByName (const MaterialPropertyName &name, MaterialData &material_data)

const MaterialProperty< T > &	getMaterialPropertyOlderByName (const MaterialPropertyName &name)

const MaterialProperty< T > &	getMaterialPropertyOlderByName (const MaterialPropertyName &name)

Moose::Kokkos::MaterialProperty< T, dimension >	getKokkosMaterialPropertyByName (const std::string &prop_name_in)

Moose::Kokkos::MaterialProperty< T, dimension >	getKokkosMaterialPropertyOldByName (const std::string &prop_name)

Moose::Kokkos::MaterialProperty< T, dimension >	getKokkosMaterialPropertyOlderByName (const std::string &prop_name)

Moose::Kokkos::MaterialProperty< T, dimension >	getKokkosMaterialProperty (const std::string &name)

Moose::Kokkos::MaterialProperty< T, dimension >	getKokkosMaterialPropertyOld (const std::string &name)

Moose::Kokkos::MaterialProperty< T, dimension >	getKokkosMaterialPropertyOlder (const std::string &name)

std::pair< const MaterialProperty< T > *, std::set< SubdomainID > >	getBlockMaterialProperty (const MaterialPropertyName &name)

const GenericMaterialProperty< T, is_ad > &	getGenericZeroMaterialProperty (const std::string &name)

const GenericMaterialProperty< T, is_ad > &	getGenericZeroMaterialProperty ()

const GenericMaterialProperty< T, is_ad > &	getGenericZeroMaterialPropertyByName (const std::string &prop_name)

const MaterialProperty< T > &	getZeroMaterialProperty (Ts... args)

std::set< SubdomainID >	getMaterialPropertyBlocks (const std::string &name)

std::vector< SubdomainName >	getMaterialPropertyBlockNames (const std::string &name)

std::set< BoundaryID >	getMaterialPropertyBoundaryIDs (const std::string &name)

std::vector< BoundaryName >	getMaterialPropertyBoundaryNames (const std::string &name)

void	checkBlockAndBoundaryCompatibility (std::shared_ptr< MaterialBase > discrete)

std::unordered_map< SubdomainID, std::vector< MaterialBase *> >	buildRequiredMaterials (bool allow_stateful=true)

void	statefulPropertiesAllowed (bool)

virtual bool	getMaterialPropertyCalled () const

virtual const std::unordered_set< unsigned int > &	getMatPropDependencies () const

virtual void	resolveOptionalProperties ()

const GenericMaterialProperty< T, is_ad > &	getPossiblyConstantGenericMaterialPropertyByName (const MaterialPropertyName &prop_name, MaterialData &material_data, const unsigned int state)

bool	isImplicit ()

Moose::StateArg	determineState () const

virtual void	threadJoin (const UserObject &) override

virtual void	threadJoin (const UserObject &) override

virtual void	subdomainSetup () override

virtual void	subdomainSetup () override

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObject (const std::string &param_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

bool	hasUserObjectByName (const UserObjectName &object_name) const

const GenericOptionalMaterialProperty< T, is_ad > &	getGenericOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const GenericOptionalMaterialProperty< T, is_ad > &	getGenericOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const OptionalMaterialProperty< T > &	getOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const OptionalMaterialProperty< T > &	getOptionalMaterialProperty (const std::string &name, const unsigned int state=0)

const OptionalADMaterialProperty< T > &	getOptionalADMaterialProperty (const std::string &name)

const OptionalADMaterialProperty< T > &	getOptionalADMaterialProperty (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOld (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOld (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOlder (const std::string &name)

const OptionalMaterialProperty< T > &	getOptionalMaterialPropertyOlder (const std::string &name)

MaterialBase &	getMaterial (const std::string &name)

MaterialBase &	getMaterial (const std::string &name)

MaterialBase &	getMaterialByName (const std::string &name, bool no_warn=false)

MaterialBase &	getMaterialByName (const std::string &name, bool no_warn=false)

bool	hasMaterialProperty (const std::string &name)

bool	hasMaterialProperty (const std::string &name)

bool	hasMaterialPropertyByName (const std::string &name)

bool	hasMaterialPropertyByName (const std::string &name)

bool	hasADMaterialProperty (const std::string &name)

bool	hasADMaterialProperty (const std::string &name)

bool	hasADMaterialPropertyByName (const std::string &name)

bool	hasADMaterialPropertyByName (const std::string &name)

bool	hasKokkosMaterialProperty (const std::string &name)

bool	hasKokkosMaterialProperty (const std::string &name)

bool	hasKokkosMaterialPropertyByName (const std::string &name)

bool	hasKokkosMaterialPropertyByName (const std::string &name)

bool	hasGenericMaterialProperty (const std::string &name)

bool	hasGenericMaterialProperty (const std::string &name)

bool	hasGenericMaterialPropertyByName (const std::string &name)

bool	hasGenericMaterialPropertyByName (const std::string &name)

const Function &	getFunction (const std::string &name) const

const Function &	getFunctionByName (const FunctionName &name) const

bool	hasFunction (const std::string &param_name) const

bool	hasFunctionByName (const FunctionName &name) const

bool	isDefaultPostprocessorValue (const std::string &param_name, const unsigned int index=0) const

bool	hasPostprocessor (const std::string &param_name, const unsigned int index=0) const

bool	hasPostprocessorByName (const PostprocessorName &name) const

std::size_t	coupledPostprocessors (const std::string &param_name) const

const PostprocessorName &	getPostprocessorName (const std::string &param_name, const unsigned int index=0) const

const VectorPostprocessorValue &	getVectorPostprocessorValue (const std::string &param_name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValue (const std::string &param_name, const std::string &vector_name, bool needs_broadcast) const

const VectorPostprocessorValue &	getVectorPostprocessorValueByName (const VectorPostprocessorName &name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValueByName (const VectorPostprocessorName &name, const std::string &vector_name, bool needs_broadcast) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOld (const std::string &param_name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOld (const std::string &param_name, const std::string &vector_name, bool needs_broadcast) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOldByName (const VectorPostprocessorName &name, const std::string &vector_name) const

const VectorPostprocessorValue &	getVectorPostprocessorValueOldByName (const VectorPostprocessorName &name, const std::string &vector_name, bool needs_broadcast) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValue (const std::string &param_name, const std::string &vector_name) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValueByName (const VectorPostprocessorName &name, const std::string &vector_name) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValueOld (const std::string &param_name, const std::string &vector_name) const

const ScatterVectorPostprocessorValue &	getScatterVectorPostprocessorValueOldByName (const VectorPostprocessorName &name, const std::string &vector_name) const

bool	hasVectorPostprocessor (const std::string &param_name, const std::string &vector_name) const

bool	hasVectorPostprocessor (const std::string &param_name) const

bool	hasVectorPostprocessorByName (const VectorPostprocessorName &name, const std::string &vector_name) const

bool	hasVectorPostprocessorByName (const VectorPostprocessorName &name) const

const VectorPostprocessorName &	getVectorPostprocessorName (const std::string &param_name) const

T &	getSampler (const std::string &name)

Sampler &	getSampler (const std::string &name)

T &	getSamplerByName (const SamplerName &name)

Sampler &	getSamplerByName (const SamplerName &name)

virtual void	meshChanged ()

virtual void	meshDisplaced ()

PerfGraph &	perfGraph ()

const PostprocessorValue &	getPostprocessorValue (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValue (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOld (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOld (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOlder (const std::string &param_name, const unsigned int index=0) const

const PostprocessorValue &	getPostprocessorValueOlder (const std::string &param_name, const unsigned int index=0) const

virtual const PostprocessorValue &	getPostprocessorValueByName (const PostprocessorName &name) const

virtual const PostprocessorValue &	getPostprocessorValueByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOldByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOldByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOlderByName (const PostprocessorName &name) const

const PostprocessorValue &	getPostprocessorValueOlderByName (const PostprocessorName &name) const

bool	isVectorPostprocessorDistributed (const std::string &param_name) const

bool	isVectorPostprocessorDistributed (const std::string &param_name) const

bool	isVectorPostprocessorDistributedByName (const VectorPostprocessorName &name) const

bool	isVectorPostprocessorDistributedByName (const VectorPostprocessorName &name) const

const Distribution &	getDistribution (const std::string &name) const

const T &	getDistribution (const std::string &name) const

const Distribution &	getDistribution (const std::string &name) const

const T &	getDistribution (const std::string &name) const

const Distribution &	getDistributionByName (const DistributionName &name) const

const T &	getDistributionByName (const std::string &name) const

const Distribution &	getDistributionByName (const DistributionName &name) const

const T &	getDistributionByName (const std::string &name) const

const Parallel::Communicator &	comm () const

processor_id_type	n_processors () const

processor_id_type	processor_id () const


void	fugacityCoefficientsHighTemp (const ADReal &pressure, const ADReal &temperature, const ADReal &co2_density, const ADReal &xco2, const ADReal &yh2o, ADReal &fco2, ADReal &fh2o) const
	Fugacity coefficients for H2O and CO2 at elevated temperatures (100C < T <= 300C). More...

ADReal	fugacityCoefficientH2OHighTemp (const ADReal &pressure, const ADReal &temperature, const ADReal &co2_density, const ADReal &xco2, const ADReal &yh2o) const

ADReal	fugacityCoefficientCO2HighTemp (const ADReal &pressure, const ADReal &temperature, const ADReal &co2_density, const ADReal &xco2, const ADReal &yh2o) const

Static Public Member Functions
static InputParameters	validParams ()

static void	callMooseError (MooseApp const app, const InputParameters &params, std::string msg, const bool with_prefix, const hit::Node node)

static void	sort (typename std::vector< T > &vector)

static void	sortDFS (typename std::vector< T > &vector)

static void	cyclicDependencyError (CyclicDependencyException< T2 > &e, const std::string &header)

Public Attributes
	usingCombinedWarningSolutionWarnings

const ConsoleStream	_console

Static Public Attributes
static const std::string	type_param

static const std::string	name_param

static const std::string	unique_name_param

static const std::string	app_param

static const std::string	moose_base_param

static const std::string	kokkos_object_param

static constexpr PropertyValue::id_type	default_property_id

static constexpr PropertyValue::id_type	zero_property_id

static constexpr auto	SYSTEM

static constexpr auto	NAME

Protected Member Functions
virtual void	checkVariables (Real pressure, Real temperature) const
	Check the input variables. More...

void	smoothCubicInterpolation (Real temperature, Real f0, Real df0, Real f1, Real df1, Real &value, Real &deriv) const
	Cubic function to smoothly interpolate between the low temperature and elevated temperature models for 99C < T < 109C. More...

virtual void	addPostprocessorDependencyHelper (const PostprocessorName &name) const override

virtual void	addVectorPostprocessorDependencyHelper (const VectorPostprocessorName &name) const override

virtual void	addUserObjectDependencyHelper (const UserObject &uo) const override

void	addReporterDependencyHelper (const ReporterName &reporter_name) override

void	flagInvalidSolutionInternal (const InvalidSolutionID invalid_solution_id) const

InvalidSolutionID	registerInvalidSolutionInternal (const std::string &message, const bool warning) const

const ReporterContextBase &	getReporterContextBaseByName (const ReporterName &reporter_name) const

const ReporterName &	getReporterName (const std::string &param_name) const

T &	declareRestartableData (const std::string &data_name, Args &&... args)

ManagedValue< T >	declareManagedRestartableDataWithContext (const std::string &data_name, void *context, Args &&... args)

const T &	getRestartableData (const std::string &data_name) const

T &	declareRestartableDataWithContext (const std::string &data_name, void *context, Args &&... args)

T &	declareRecoverableData (const std::string &data_name, Args &&... args)

T &	declareRestartableDataWithObjectName (const std::string &data_name, const std::string &object_name, Args &&... args)

T &	declareRestartableDataWithObjectNameWithContext (const std::string &data_name, const std::string &object_name, void *context, Args &&... args)

std::string	restartableName (const std::string &data_name) const

const T &	getMeshProperty (const std::string &data_name, const std::string &prefix)

const T &	getMeshProperty (const std::string &data_name)

bool	hasMeshProperty (const std::string &data_name, const std::string &prefix) const

bool	hasMeshProperty (const std::string &data_name, const std::string &prefix) const

bool	hasMeshProperty (const std::string &data_name) const

bool	hasMeshProperty (const std::string &data_name) const

std::string	meshPropertyName (const std::string &data_name) const

PerfID	registerTimedSection (const std::string &section_name, const unsigned int level) const

PerfID	registerTimedSection (const std::string &section_name, const unsigned int level, const std::string &live_message, const bool print_dots=true) const

std::string	timedSectionName (const std::string &section_name) const

bool	isCoupledScalar (const std::string &var_name, unsigned int i=0) const

unsigned int	coupledScalarComponents (const std::string &var_name) const

unsigned int	coupledScalar (const std::string &var_name, unsigned int comp=0) const

libMesh::Order	coupledScalarOrder (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarValue (const std::string &var_name, unsigned int comp=0) const

const ADVariableValue &	adCoupledScalarValue (const std::string &var_name, unsigned int comp=0) const

const GenericVariableValue< is_ad > &	coupledGenericScalarValue (const std::string &var_name, unsigned int comp=0) const

const GenericVariableValue< false > &	coupledGenericScalarValue (const std::string &var_name, const unsigned int comp) const

const GenericVariableValue< true > &	coupledGenericScalarValue (const std::string &var_name, const unsigned int comp) const

const VariableValue &	coupledVectorTagScalarValue (const std::string &var_name, TagID tag, unsigned int comp=0) const

const VariableValue &	coupledMatrixTagScalarValue (const std::string &var_name, TagID tag, unsigned int comp=0) const

const VariableValue &	coupledScalarValueOld (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarValueOlder (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDot (const std::string &var_name, unsigned int comp=0) const

const ADVariableValue &	adCoupledScalarDot (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDot (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotOld (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDotOld (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDu (const std::string &var_name, unsigned int comp=0) const

const VariableValue &	coupledScalarDotDotDu (const std::string &var_name, unsigned int comp=0) const

const MooseVariableScalar *	getScalarVar (const std::string &var_name, unsigned int comp) const

virtual void	checkMaterialProperty (const std::string &name, const unsigned int state)

virtual void	getKokkosMaterialPropertyHook (const std::string &, const unsigned int)

void	markMatPropRequested (const std::string &)

MaterialPropertyName	getMaterialPropertyName (const std::string &name) const

void	checkExecutionStage ()

const T &	getReporterValue (const std::string &param_name, const std::size_t time_index=0)

const T &	getReporterValue (const std::string &param_name, ReporterMode mode, const std::size_t time_index=0)

const T &	getReporterValue (const std::string &param_name, const std::size_t time_index=0)

const T &	getReporterValue (const std::string &param_name, ReporterMode mode, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, ReporterMode mode, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, const std::size_t time_index=0)

const T &	getReporterValueByName (const ReporterName &reporter_name, ReporterMode mode, const std::size_t time_index=0)

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValue (const std::string &param_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

bool	hasReporterValueByName (const ReporterName &reporter_name) const

const GenericMaterialProperty< T, is_ad > *	defaultGenericMaterialProperty (const std::string &name)

const GenericMaterialProperty< T, is_ad > *	defaultGenericMaterialProperty (const std::string &name)

const MaterialProperty< T > *	defaultMaterialProperty (const std::string &name)

const MaterialProperty< T > *	defaultMaterialProperty (const std::string &name)

const ADMaterialProperty< T > *	defaultADMaterialProperty (const std::string &name)

const ADMaterialProperty< T > *	defaultADMaterialProperty (const std::string &name)


void	funcABHighTemp (Real pressure, Real temperature, Real Xnacl, Real co2_density, Real xco2, Real yh2o, Real &A, Real &B) const
	The function A (Eq. More...

void	funcABHighTemp (Real pressure, Real temperature, Real Xnacl, Real co2_density, Real xco2, Real yh2o, Real &A, Real &dA_dp, Real &dA_dT, Real &B, Real &dB_dp, Real &dB_dT, Real &dB_dX) const

void	funcABLowTemp (const ADReal &pressure, const ADReal &temperature, const ADReal &co2_density, ADReal &A, ADReal &B) const

Static Protected Member Functions
static std::string	meshPropertyName (const std::string &data_name, const std::string &prefix)

Protected Attributes
const unsigned int	_salt_component
	Salt component index. More...

const BrineFluidProperties &	_brine_fp
	Fluid properties UserObject for water. More...

const SinglePhaseFluidProperties &	_co2_fp
	Fluid properties UserObject for the CO2. More...

const Real	_Mh2o
	Molar mass of water (kg/mol) More...

const Real	_invMh2o
	Inverse of molar mass of H2O (mol/kg) More...

const Real	_Mco2
	Molar mass of CO2 (kg/mol) More...

const Real	_Mnacl
	Molar mass of NaCL. More...

const Real	_Rbar
	Molar gas constant in bar cm^3 /(K mol) More...

const Real	_Tlower
	Temperature below which the Spycher, Pruess & Ennis-King (2003) model is used (K) More...

const Real	_Tupper
	Temperature above which the Spycher & Pruess (2010) model is used (K) More...

const Real	_Zmin
	Minimum Z - below this value all CO2 will be dissolved. More...

const std::vector< Real >	_co2_henry
	Henry's coefficeients for CO2. More...

const unsigned int	_aqueous_fluid_component
	Fluid component number of the aqueous component. More...

unsigned int	_gas_fluid_component
	Fluid component number of the gas phase. More...

const unsigned int	_pidx
	Index of derivative wrt pressure. More...

const unsigned int	_Zidx
	Index of derivative wrt total mass fraction Z. More...

const unsigned int	_Tidx
	Index of derivative wrt temperature. More...

const unsigned int	_Xidx
	Index of derivative wrt salt mass fraction X. More...

const Real	_nr_max_its
	Maximum number of iterations for the Newton-Raphson routine. More...

const Real	_nr_tol
	Tolerance for Newton-Raphson iterations. More...

unsigned int	_num_phases
	Number of phases. More...

unsigned int	_num_components
	Number of components. More...

const unsigned int	_aqueous_phase_number
	Phase number of the aqueous phase. More...

unsigned int	_gas_phase_number
	Phase number of the gas phase. More...

const Real	_R
	Universal gas constant (J/mol/K) More...

const Real	_T_c2k
	Conversion from C to K. More...

const PorousFlowCapillaryPressure &	_pc
	Capillary pressure UserObject. More...

FluidStateProperties	_empty_fsp
	Empty FluidStateProperties object. More...

SubProblem &	_subproblem

FEProblemBase &	_fe_problem

SystemBase &	_sys

const THREAD_ID	_tid

Assembly &	_assembly

const Moose::CoordinateSystemType &	_coord_sys

const bool	_duplicate_initial_execution

std::set< std::string >	_depend_uo

const bool &	_enabled

MooseApp &	_app

Factory &	_factory

ActionFactory &	_action_factory

const std::string &	_type

const std::string &	_name

const InputParameters &	_pars

const ExecFlagEnum &	_execute_enum

const ExecFlagType &	_current_execute_flag

MooseApp &	_restartable_app

const std::string	_restartable_system_name

const THREAD_ID	_restartable_tid

const bool	_restartable_read_only

FEProblemBase &	_mci_feproblem

FEProblemBase &	_mdi_feproblem

MooseApp &	_pg_moose_app

const std::string	_prefix

FEProblemBase &	_sc_fe_problem

const THREAD_ID	_sc_tid

const Real &	_real_zero

const VariableValue &	_scalar_zero

const Point &	_point_zero

const InputParameters &	_mi_params

const std::string	_mi_name

const MooseObjectName	_mi_moose_object_name

FEProblemBase &	_mi_feproblem

SubProblem &	_mi_subproblem

const THREAD_ID	_mi_tid

const bool	_is_kokkos_object

const Moose::MaterialDataType	_material_data_type

MaterialData &	_material_data

bool	_stateful_allowed

bool	_get_material_property_called

std::vector< std::unique_ptr< PropertyValue > >	_default_properties

std::unordered_set< unsigned int >	_material_property_dependencies

const MaterialPropertyName	_get_suffix

const bool	_use_interpolated_state

const InputParameters &	_ti_params

FEProblemBase &	_ti_feproblem

bool	_is_implicit

Real &	_t

const Real &	_t_old

int &	_t_step

Real &	_dt

Real &	_dt_old

bool	_is_transient

const Parallel::Communicator &	_communicator

Static Protected Attributes
static const std::string	_interpolated_old

static const std::string	_interpolated_older

Detailed Description

Specialized class for brine and CO2 including calculation of mutual solubility of the two fluids using a high-accuracy fugacity-based formulation.

For temperatures 12C <= T <= 99C, the formulation is based on Spycher, Pruess and Ennis-King, CO2-H2O mixtures in the geological sequestration of CO2. I. Assessment and calculation of mutual solubilities from 12 to 100C and up to 600 bar, Geochimica et Cosmochimica Acta, 67, 3015-3031 (2003), and Spycher and Pruess, CO2-H2O mixtures in the geological sequestration of CO2. II. Partitioning in chloride brine at 12-100C and up to 600 bar, Geochimica et Cosmochimica Acta, 69, 3309-3320 (2005).

For temperatures 109C <= T <= 300C, the formulation is based on Spycher and Pruess, A Phase-Partitioning Model for CO2-Brine Mixtures at Elevated Temperatures and Pressures: Application to CO2-Enhanced Geothermal Systems, Transport in Porous Media, 82, 173-196 (2010)

As the two formulations do not coincide at temperatures near 100C, a cubic polynomial is used in the intermediate temperature range 99C < T < 109C to provide a smooth transition from the two formulations in this region.

Notation convention Throughout this class, both mole fractions and mass fractions will be used. The following notation will be used: yk: mole fraction of component k in the gas phase xk: mole fraction of component k in the liquid phase Yk: mass fraction of component k in the gas phase Xk: mass fraction of component k in the liquid phase

Definition at line 47 of file PorousFlowBrineCO2.h.

Constructor & Destructor Documentation

◆ PorousFlowBrineCO2()

PorousFlowBrineCO2::PorousFlowBrineCO2 ( const InputParameters & parameters )

Definition at line 29 of file PorousFlowBrineCO2.C.

   : PorousFlowFluidStateMultiComponentBase(parameters),
     _salt_component(getParam<unsigned int>("salt_component")),
     _brine_fp(getUserObject<BrineFluidProperties>("brine_fp")),
     _co2_fp(getUserObject<SinglePhaseFluidProperties>("co2_fp")),
     _Mh2o(_brine_fp.molarMassH2O()),
     _invMh2o(1.0 / _Mh2o),
     _Mco2(_co2_fp.molarMass()),
     _Mnacl(_brine_fp.molarMassNaCl()),
     _Rbar(_R * 10.0),
     _Tlower(372.15),
     _Tupper(382.15),
     _Zmin(1.0e-4),
     _co2_henry(_co2_fp.henryCoefficients())
 {
   // Check that the correct FluidProperties UserObjects have been provided
   if (_co2_fp.fluidName() != "co2")
     paramError("co2_fp", "A valid CO2 FluidProperties UserObject must be provided");
 
   if (_brine_fp.fluidName() != "brine")
     paramError("brine_fp", "A valid Brine FluidProperties UserObject must be provided");
 
   // Set the number of phases and components, and their indexes
   _num_phases = 2;
   _num_components = 3;
   _gas_phase_number = 1 - _aqueous_phase_number;
   _gas_fluid_component = 3 - _aqueous_fluid_component - _salt_component;
 
   // Check that _aqueous_phase_number is <= total number of phases
   if (_aqueous_phase_number >= _num_phases)
     paramError("liquid_phase_number",
                "This value is larger than the possible number of phases ",
                _num_phases);
 
   // Check that _aqueous_fluid_component is <= total number of fluid components
   if (_aqueous_fluid_component >= _num_components)
     paramError("liquid_fluid_component",
                "This value is larger than the possible number of fluid components",
                _num_components);
 
   // Check that the salt component index is not identical to the liquid fluid component
   if (_salt_component == _aqueous_fluid_component)
     paramError(
         "salt_component",
         "The value provided must be different from the value entered in liquid_fluid_component");
 
   // Check that _salt_component is <= total number of fluid components
   if (_salt_component >= _num_components)
     paramError("salt_component",
                "The value provided is larger than the possible number of fluid components",
                _num_components);
 
   _empty_fsp = FluidStateProperties(_num_components);
 }

Member Function Documentation

◆ activityCoefficient()

ADReal PorousFlowBrineCO2::activityCoefficient	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl
	)		const

Activity coefficient for CO2 in brine.

From Duan and Sun, An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 257 to 533 K and from 0 to 2000 bar, Chem. Geol., 193, 257-271 (2003)

Note: this is not a 'true' activity coefficient, and is instead related to the molality of CO2 in water and brine. Nevertheless, Spycher and Pruess (2005) refer to it as an activity coefficient, so this notation is followed here.

Parameters

pressure	phase pressure (Pa)
temperature	phase temperature (K)
Xnacl	salt mass fraction (kg/kg)

Returns: activity coefficient for CO2 in brine

Definition at line 588 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMoleFractionsLowTemp().

 {
   using std::log, std::exp;
 
   // Need pressure in bar
   const ADReal pbar = pressure * 1.0e-5;
   // Need NaCl molality (mol/kg)
   const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
 
   const ADReal lambda = -0.411370585 + 6.07632013e-4 * temperature + 97.5347708 / temperature -
                         0.0237622469 * pbar / temperature +
                         0.0170656236 * pbar / (630.0 - temperature) +
                         1.41335834e-5 * temperature * log(pbar);
 
   const ADReal xi = 3.36389723e-4 - 1.9829898e-5 * temperature +
                     2.12220830e-3 * pbar / temperature -
                     5.24873303e-3 * pbar / (630.0 - temperature);
 
   return exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
 }

◆ activityCoefficientCO2()

ADReal PorousFlowBrineCO2::activityCoefficientCO2	(	const ADReal &	temperature,
		const ADReal &	xco2
	)		const

Activity coefficient of CO2 Eq.

(13) from Spycher & Pruess (2010)

Parameters

temperature	temperature (K)
xco2	mole fraction of CO2 in liquid phase (-)

Returns: activity coefficient

Definition at line 571 of file PorousFlowBrineCO2.C.

Referenced by funcABHighTemp().

 {
   using std::exp;
 
   if (temperature.value() <= 373.15)
     return 1.0;
   else
   {
     const ADReal Tref = temperature - 373.15;
     const ADReal xh2o = 1.0 - xco2;
     const ADReal Am = -3.084e-2 * Tref + 1.927e-5 * Tref * Tref;
 
     return exp(2.0 * Am * xco2 * xh2o * xh2o);
   }
 }

◆ activityCoefficientH2O()

ADReal PorousFlowBrineCO2::activityCoefficientH2O	(	const ADReal &	temperature,
		const ADReal &	xco2
	)		const

Activity coefficient of H2O Eq.

(12) from Spycher & Pruess (2010)

Parameters

temperature	temperature (K)
xco2	mole fraction of CO2 in liquid phase (-)

Returns: activity coefficient

Definition at line 554 of file PorousFlowBrineCO2.C.

Referenced by funcABHighTemp().

 {
   using std::exp;
 
   if (temperature.value() <= 373.15)
     return 1.0;
   else
   {
     const ADReal Tref = temperature - 373.15;
     const ADReal xh2o = 1.0 - xco2;
     const ADReal Am = -3.084e-2 * Tref + 1.927e-5 * Tref * Tref;
 
     return exp((Am - 2.0 * Am * xh2o) * xco2 * xco2);
   }
 }

◆ activityCoefficientHighTemp()

ADReal PorousFlowBrineCO2::activityCoefficientHighTemp	(	const ADReal &	temperature,
		const ADReal &	Xnacl
	)		const

Activity coefficient for CO2 in brine used in the elevated temperature formulation.

Eq. (18) from Spycher and Pruess (2010).

Note: unlike activityCoefficient(), no pressure dependence is included in this formulation

Parameters

temperature	phase temperature (K)
Xnacl	salt mass fraction (kg/kg)

Returns: activity coefficient for CO2 in brine

Definition at line 612 of file PorousFlowBrineCO2.C.

Referenced by funcABHighTemp().

 {
   using std::exp;
 
   // Need NaCl molality (mol/kg)
   const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
 
   const ADReal T2 = temperature * temperature;
   const ADReal T3 = temperature * T2;
 
   const ADReal lambda = 2.217e-4 * temperature + 1.074 / temperature + 2648.0 / T2;
   const ADReal xi = 1.3e-5 * temperature - 20.12 / temperature + 5259.0 / T2;
 
   return (1.0 - mnacl / _invMh2o) * exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
 }

◆ aqueousComponentIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::aqueousComponentIndex ( ) const

inlineinherited

The index of the aqueous fluid component.

Returns: aqueous fluid component number

Definition at line 29 of file PorousFlowFluidStateMultiComponentBase.h.

29 { return _aqueous_fluid_component; };

PorousFlowFluidStateMultiComponentBase::_aqueous_fluid_component

const unsigned int _aqueous_fluid_component

Fluid component number of the aqueous component.

Definition: PorousFlowFluidStateMultiComponentBase.h:95

◆ aqueousPhaseIndex()

unsigned int PorousFlowFluidStateBase::aqueousPhaseIndex ( ) const

inlineinherited

The index of the aqueous phase.

Returns: aqueous phase number

Definition at line 77 of file PorousFlowFluidStateBase.h.

77 { return _aqueous_phase_number; };

PorousFlowFluidStateBase::_aqueous_phase_number

const unsigned int _aqueous_phase_number

Phase number of the aqueous phase.

Definition: PorousFlowFluidStateBase.h:102

◆ checkVariables()

void PorousFlowBrineCO2::checkVariables	(	Real	pressure,
		Real	temperature
	)		const

protectedvirtual

Check the input variables.

Definition at line 1181 of file PorousFlowBrineCO2.C.

Referenced by thermophysicalProperties(), and totalMassFraction().

 {
   // The calculation of mass fractions is valid from 12C <= T <= 300C, and
   // pressure less than 60 MPa
   if (temperature < 285.15 || temperature > 573.15)
     mooseException(name() + ": temperature " + Moose::stringify(temperature) +
                    " is outside range 285.15 K <= T <= 573.15 K");
 
   if (pressure > 6.0e7)
     mooseException(name() + ": pressure " + Moose::stringify(pressure) +
                    " must be less than 60 MPa");
 }

◆ clearFluidStateProperties()

void PorousFlowFluidStateBase::clearFluidStateProperties ( std::vector< FluidStateProperties > & fsp ) const

inherited

Clears the contents of the FluidStateProperties data structure.

Parameters

[out] fsp FluidStateProperties data structure with all data initialized to 0

Definition at line 33 of file PorousFlowFluidStateBase.C.

Referenced by PorousFlowWaterNCG::thermophysicalProperties(), and thermophysicalProperties().

 {
   std::fill(fsp.begin(), fsp.end(), _empty_fsp);
 }

◆ enthalpyOfDissolution()

ADReal PorousFlowBrineCO2::enthalpyOfDissolution ( const ADReal & temperature ) const

Enthalpy of dissolution of CO2 in brine calculated using linear fit to model of Duan and Sun, An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar, Chemical geology, 193, 257–271 (2003).

In the region of interest, the more complicated model given in Eq. (8) of Duan and Sun is well approximated by a simple linear fit (R^2 = 0.9997).

Note: as the effect of salt mass fraction is small, it is not included in this function.

Parameters

temperature fluid temperature (K)

Returns: enthalpy of dissolution (J/kg)

Definition at line 1152 of file PorousFlowBrineCO2.C.

Referenced by liquidProperties().

 {
   // Linear fit to model of Duan and Sun (2003) (in kJ/mol)
   const ADReal delta_h = -58.3533 + 0.134519 * temperature;
 
   // Convert to J/kg
   return delta_h * 1000.0 / _Mco2;
 }

◆ enthalpyOfDissolutionGas()

ADReal PorousFlowBrineCO2::enthalpyOfDissolutionGas	(	const ADReal &	temperature,
		const ADReal &	Xnacl
	)		const

Enthalpy of dissolution of gas phase CO2 in brine calculated using Henry's constant From Himmelblau, Partial molal heats and entropies of solution for gases dissolved in water from the freezing to the near critical point, J.

Phys. Chem. 63 (1959). Correction due to salinity from Battistelli et al, A fluid property module for the TOUGH2 simulator for saline brines with non-condensible gas, Proc. Eighteenth Workshop on Geothermal Reservoir Engineering (1993).

Parameters

temperature	fluid temperature (K)
Xnacl	NaCl mass fraction (kg/kg)

Returns: enthalpy of dissolution (J/kg)

Definition at line 1122 of file PorousFlowBrineCO2.C.

 {
   // Henry's constant
   const ADReal Kh = henryConstant(temperature, Xnacl);
 
   ADReal hdis = -_R * temperature * temperature * Kh.derivatives()[_Tidx] / Kh / _Mco2;
 
   // Derivative of enthalpy of dissolution wrt temperature and xnacl requires the second
   // derivatives of Henry's constant. For simplicity, approximate these numerically
   const Real dT = temperature.value() * 1.0e-8;
   const ADReal T2 = temperature + dT;
   const ADReal Kh2 = henryConstant(T2, Xnacl);
 
   const Real dhdis_dT =
       (-_R * T2 * T2 * Kh2.derivatives()[_Tidx] / Kh2 / _Mco2 - hdis).value() / dT;
 
   const Real dX = Xnacl.value() * 1.0e-8;
   const ADReal X2 = Xnacl + dX;
   const ADReal Kh3 = henryConstant(temperature, X2);
 
   const Real dhdis_dX =
       (-_R * temperature * temperature * Kh3.derivatives()[_Tidx] / Kh3 / _Mco2 - hdis).value() /
       dX;
 
   hdis.derivatives() = temperature.derivatives() * dhdis_dT + Xnacl.derivatives() * dhdis_dX;
 
   return hdis;
 }

◆ equilibriumConstantCO2()

ADReal PorousFlowBrineCO2::equilibriumConstantCO2 ( const ADReal & temperature ) const

Equilibrium constant for CO2 For temperatures 12C <= T <= 99C, uses Spycher, Pruess and Ennis-King (2003) For temperatures 109C <= T <= 300C, uses Spycher and Pruess (2010) For temperatures 99C < T < 109C, the value is calculated by smoothly interpolating the two formulations.

Parameters

temperature temperature (K)

Returns: equilibrium constant for CO2

Definition at line 659 of file PorousFlowBrineCO2.C.

Referenced by funcABHighTemp(), and funcABLowTemp().

 {
   using std::pow;
 
   // Uses temperature in Celsius
   const ADReal Tc = temperature - _T_c2k;
   const ADReal Tc2 = Tc * Tc;
   const ADReal Tc3 = Tc2 * Tc;
 
   ADReal logK0CO2;
 
   if (Tc <= 99.0)
     logK0CO2 = 1.189 + 1.304e-2 * Tc - 5.446e-5 * Tc2;
 
   else if (Tc > 99.0 && Tc < 109.0)
   {
     const ADReal Tint = (Tc - 99.0) / 10.0;
     const ADReal Tint2 = Tint * Tint;
     logK0CO2 = 1.9462 + 2.25692e-2 * Tint - 9.49577e-3 * Tint2 - 6.77721e-3 * Tint * Tint2;
   }
 
   else // 109 <= Tc <= 300
     logK0CO2 = 1.668 + 3.992e-3 * Tc - 1.156e-5 * Tc2 + 1.593e-9 * Tc3;
 
   return pow(10.0, logK0CO2);
 }

◆ equilibriumConstantH2O()

ADReal PorousFlowBrineCO2::equilibriumConstantH2O ( const ADReal & temperature ) const

Equilibrium constant for H2O For temperatures 12C <= T <= 99C, uses Spycher, Pruess and Ennis-King (2003) For temperatures 109C <= T <= 300C, uses Spycher and Pruess (2010) For temperatures 99C < T < 109C, the value is calculated by smoothly interpolating the two formulations.

Parameters

temperature temperature (K)

Returns: equilibrium constant for H2O

Definition at line 630 of file PorousFlowBrineCO2.C.

Referenced by funcABHighTemp(), and funcABLowTemp().

 {
   using std::pow;
 
   // Uses temperature in Celsius
   const ADReal Tc = temperature - _T_c2k;
   const ADReal Tc2 = Tc * Tc;
   const ADReal Tc3 = Tc2 * Tc;
   const ADReal Tc4 = Tc3 * Tc;
 
   ADReal logK0H2O;
 
   if (Tc <= 99.0)
     logK0H2O = -2.209 + 3.097e-2 * Tc - 1.098e-4 * Tc2 + 2.048e-7 * Tc3;
 
   else if (Tc > 99.0 && Tc < 109.0)
   {
     const ADReal Tint = (Tc - 99.0) / 10.0;
     const ADReal Tint2 = Tint * Tint;
     logK0H2O = -0.0204026 + 0.0152513 * Tint + 0.417565 * Tint2 - 0.278636 * Tint * Tint2;
   }
 
   else // 109 <= Tc <= 300
     logK0H2O = -2.1077 + 2.8127e-2 * Tc - 8.4298e-5 * Tc2 + 1.4969e-7 * Tc3 - 1.1812e-10 * Tc4;
 
   return pow(10.0, logK0H2O);
 }

◆ equilibriumMassFractions()

void PorousFlowBrineCO2::equilibriumMassFractions	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		ADReal &	Xco2,
		ADReal &	Yh2o
	)		const

Mass fractions of CO2 in brine and water vapor in CO2 at equilibrium.

Parameters

	pressure	phase pressure (Pa)
	temperature	phase temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
[out]	Xco2	mass fraction of CO2 in liquid (kg/kg)
[out]	Yh2o	mass fraction of H2O in gas (kg/kg)

Definition at line 372 of file PorousFlowBrineCO2.C.

Referenced by massFractions(), and totalMassFraction().

 {
   // Mole fractions at equilibrium
   ADReal xco2, yh2o;
   equilibriumMoleFractions(pressure, temperature, Xnacl, xco2, yh2o);
 
   // The mass fraction of H2O in gas (assume no salt in gas phase) and derivatives
   // wrt p, T, and X
   Yh2o = yh2o * _Mh2o / (yh2o * _Mh2o + (1.0 - yh2o) * _Mco2);
 
   // NaCl molality (mol/kg)
   const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
 
   // The molality of CO2 in 1kg of H2O
   const ADReal mco2 = xco2 * (2.0 * mnacl + _invMh2o) / (1.0 - xco2);
   // The mass fraction of CO2 in brine is then
   const ADReal denominator = (1.0 + mnacl * _Mnacl + mco2 * _Mco2);
   Xco2 = mco2 * _Mco2 / denominator;
 }

◆ equilibriumMoleFractions()

void PorousFlowBrineCO2::equilibriumMoleFractions	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		ADReal &	xco2,
		ADReal &	yh2o
	)		const

Mole fractions of CO2 in brine and water vapor in CO2 at equilibrium.

In the low temperature regime (T <= 99C), the mole fractions are calculated directly, while in the elevated temperature regime (T >= 109C), they are calculated iteratively. In the intermediate regime (99C < T < 109C), a cubic polynomial is used to smoothly connect the low and elevated temperature regimes.

Parameters

	pressure	phase pressure (Pa)
	temperature	phase temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
[out]	xco2	mole fraction of CO2 in liquid
[out]	yh2o	mole fraction of H2O in gas

Definition at line 687 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMassFractions().

 {
   if (temperature.value() <= _Tlower)
   {
     equilibriumMoleFractionsLowTemp(pressure, temperature, Xnacl, xco2, yh2o);
   }
   else if (temperature.value() > _Tlower && temperature.value() < _Tupper)
   {
     // Cubic polynomial in this regime
     const Real Tint = (temperature.value() - _Tlower) / 10.0;
 
     // Equilibrium mole fractions and derivatives at the lower temperature
     ADReal Tlower = _Tlower;
     Moose::derivInsert(Tlower.derivatives(), _Tidx, 1.0);
 
     ADReal xco2_lower, yh2o_lower;
     equilibriumMoleFractionsLowTemp(pressure, Tlower, Xnacl, xco2_lower, yh2o_lower);
 
     const Real dxco2_dT_lower = xco2_lower.derivatives()[_Tidx];
     const Real dyh2o_dT_lower = yh2o_lower.derivatives()[_Tidx];
 
     // Equilibrium mole fractions and derivatives at the upper temperature
     Real xco2_upper, yh2o_upper;
     Real co2_density_upper = _co2_fp.rho_from_p_T(pressure.value(), _Tupper);
 
     solveEquilibriumMoleFractionHighTemp(
         pressure.value(), _Tupper, Xnacl.value(), co2_density_upper, xco2_upper, yh2o_upper);
 
     Real A, dA_dp, dA_dT, B, dB_dp, dB_dT, dB_dX;
     funcABHighTemp(pressure.value(),
                    _Tupper,
                    Xnacl.value(),
                    co2_density_upper,
                    xco2_upper,
                    yh2o_upper,
                    A,
                    dA_dp,
                    dA_dT,
                    B,
                    dB_dp,
                    dB_dT,
                    dB_dX);
 
     const Real dyh2o_dT_upper =
         ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
     const Real dxco2_dT_upper = dB_dT * (1.0 - yh2o_upper) - B * dyh2o_dT_upper;
 
     // The mole fractions in this regime are then found by interpolation
     Real xco2r, yh2or, dxco2_dT, dyh2o_dT;
     smoothCubicInterpolation(
         Tint, xco2_lower.value(), dxco2_dT_lower, xco2_upper, dxco2_dT_upper, xco2r, dxco2_dT);
     smoothCubicInterpolation(
         Tint, yh2o_lower.value(), dyh2o_dT_lower, yh2o_upper, dyh2o_dT_upper, yh2or, dyh2o_dT);
 
     xco2 = xco2r;
     Moose::derivInsert(xco2.derivatives(), _pidx, xco2_lower.derivatives()[_pidx]);
     Moose::derivInsert(xco2.derivatives(), _Tidx, dxco2_dT);
     Moose::derivInsert(xco2.derivatives(), _Xidx, xco2_lower.derivatives()[_Xidx]);
 
     yh2o = yh2or;
     Moose::derivInsert(yh2o.derivatives(), _pidx, yh2o_lower.derivatives()[_pidx]);
     Moose::derivInsert(yh2o.derivatives(), _Tidx, dyh2o_dT);
     Moose::derivInsert(yh2o.derivatives(), _Xidx, yh2o_lower.derivatives()[_Xidx]);
   }
   else
   {
     // CO2 density and derivatives wrt pressure and temperature
     const Real co2_density = _co2_fp.rho_from_p_T(pressure.value(), temperature.value());
 
     // Equilibrium mole fractions solved using iteration in this regime
     Real xco2r, yh2or;
     solveEquilibriumMoleFractionHighTemp(
         pressure.value(), temperature.value(), Xnacl.value(), co2_density, xco2r, yh2or);
 
     // Can use these in funcABHighTemp() to get derivatives analytically rather than by iteration
     Real A, dA_dp, dA_dT, B, dB_dp, dB_dT, dB_dX;
     funcABHighTemp(pressure.value(),
                    temperature.value(),
                    Xnacl.value(),
                    co2_density,
                    xco2r,
                    yh2or,
                    A,
                    dA_dp,
                    dA_dT,
                    B,
                    dB_dp,
                    dB_dT,
                    dB_dX);
 
     const Real dyh2o_dp =
         ((1.0 - B) * dA_dp + (A - 1.0) * A * dB_dp) / (1.0 - A * B) / (1.0 - A * B);
     const Real dxco2_dp = dB_dp * (1.0 - yh2or) - B * dyh2o_dp;
 
     const Real dyh2o_dT =
         ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
     const Real dxco2_dT = dB_dT * (1.0 - yh2or) - B * dyh2o_dT;
 
     const Real dyh2o_dX = ((A - 1.0) * A * dB_dX) / (1.0 - A * B) / (1.0 - A * B);
     const Real dxco2_dX = dB_dX * (1.0 - yh2or) - B * dyh2o_dX;
 
     xco2 = xco2r;
     Moose::derivInsert(xco2.derivatives(), _pidx, dxco2_dp);
     Moose::derivInsert(xco2.derivatives(), _Tidx, dxco2_dT);
     Moose::derivInsert(xco2.derivatives(), _Xidx, dxco2_dX);
 
     yh2o = yh2or;
     Moose::derivInsert(yh2o.derivatives(), _pidx, dyh2o_dp);
     Moose::derivInsert(yh2o.derivatives(), _Tidx, dyh2o_dT);
     Moose::derivInsert(yh2o.derivatives(), _Xidx, dyh2o_dX);
   }
 }

◆ equilibriumMoleFractionsLowTemp()

void PorousFlowBrineCO2::equilibriumMoleFractionsLowTemp	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		ADReal &	xco2,
		ADReal &	yh2o
	)		const

Mole fractions of CO2 in brine and water vapor in CO2 at equilibrium in the low temperature regime (T <= 99C).

Parameters

	pressure	phase pressure (Pa)
	temperature	phase temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
[out]	xco2	mole fraction of CO2 in liquid
[out]	yh2o	mass fraction of mole in gas

Definition at line 805 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMoleFractions().

 {
   if (temperature.value() > 373.15)
     mooseError(name(),
                ": equilibriumMoleFractionsLowTemp() is not valid for T > 373.15K. Use "
                "equilibriumMoleFractions() instead");
 
   // CO2 density and derivatives wrt pressure and temperature
   const ADReal co2_density = _co2_fp.rho_from_p_T(pressure, temperature);
 
   // Assume infinite dilution (yh20 = 0 and xco2 = 0) in low temperature regime
   ADReal A, B;
   funcABLowTemp(pressure, temperature, co2_density, A, B);
 
   // As the activity coefficient for CO2 in brine used in this regime isn't a 'true'
   // activity coefficient, we instead calculate the molality of CO2 in water, then
   // correct it for brine, and then calculate the mole fractions.
   // The mole fraction in pure water is
   const ADReal yh2ow = (1.0 - B) / (1.0 / A - B);
   const ADReal xco2w = B * (1.0 - yh2ow);
 
   // Molality of CO2 in pure water
   const ADReal mco2w = xco2w * _invMh2o / (1.0 - xco2w);
   // Molality of CO2 in brine is then calculated using gamma
   const ADReal gamma = activityCoefficient(pressure, temperature, Xnacl);
   const ADReal mco2 = mco2w / gamma;
 
   // Need NaCl molality (mol/kg)
   const ADReal mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
 
   // Mole fractions of CO2 and H2O in liquid and gas phases
   const ADReal total_moles = 2.0 * mnacl + _invMh2o + mco2;
   xco2 = mco2 / total_moles;
   yh2o = A * (1.0 - xco2 - 2.0 * mnacl / total_moles);
 }

◆ execute()

void PorousFlowFluidStateBase::execute ( )

inlinefinalvirtualinherited

Implements GeneralUserObject.

Definition at line 58 of file PorousFlowFluidStateBase.h.

58 {};

◆ finalize()

void PorousFlowFluidStateBase::finalize ( )

inlinefinalvirtualinherited

Implements GeneralUserObject.

Definition at line 59 of file PorousFlowFluidStateBase.h.

59 {};

◆ fluidStateName()

std::string PorousFlowBrineCO2::fluidStateName ( ) const

overridevirtual

Name of FluidState.

Implements PorousFlowFluidStateBase.

Definition at line 85 of file PorousFlowBrineCO2.C.

 {
   return "brine-co2";
 }

◆ fugacityCoefficientCO2HighTemp()

ADReal PorousFlowBrineCO2::fugacityCoefficientCO2HighTemp	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	co2_density,
		const ADReal &	xco2,
		const ADReal &	yh2o
	)		const

Definition at line 506 of file PorousFlowBrineCO2.C.

Referenced by fugacityCoefficientsHighTemp().

 {
   using std::sqrt, std::pow, std::log, std::exp;
 
   // Need pressure in bar
   const ADReal pbar = pressure * 1.0e-5;
   // Molar volume in cm^3/mol
   const ADReal V = _Mco2 / co2_density * 1.0e6;
 
   // Redlich-Kwong parameters
   const ADReal yco2 = 1.0 - yh2o;
   const ADReal xh2o = 1.0 - xco2;
 
   const ADReal aCO2 = 8.008e7 - 4.984e4 * temperature;
   const ADReal aH2O = 1.337e8 - 1.4e4 * temperature;
   const Real bCO2 = 28.25;
   const Real bH2O = 15.7;
   const ADReal KH2OCO2 = 1.427e-2 - 4.037e-4 * temperature;
   const ADReal KCO2H2O = 0.4228 - 7.422e-4 * temperature;
   const ADReal kH2OCO2 = KH2OCO2 * yh2o + KCO2H2O * yco2;
   const ADReal kCO2H2O = KCO2H2O * yh2o + KH2OCO2 * yco2;
 
   const ADReal aH2OCO2 = sqrt(aCO2 * aH2O) * (1.0 - kH2OCO2);
   const ADReal aCO2H2O = sqrt(aCO2 * aH2O) * (1.0 - kCO2H2O);
 
   const ADReal amix = yh2o * yh2o * aH2O + yh2o * yco2 * (aH2OCO2 + aCO2H2O) + yco2 * yco2 * aCO2;
   const ADReal bmix = yh2o * bH2O + yco2 * bCO2;
 
   const ADReal t15 = pow(temperature, 1.5);
 
   ADReal lnPhiCO2 = bCO2 / bmix * (pbar * V / (_Rbar * temperature) - 1.0) -
                     log(pbar * (V - bmix) / (_Rbar * temperature));
 
   ADReal term3 = (2.0 * yco2 * aCO2 + yh2o * (aH2OCO2 + aCO2H2O) -
                   yh2o * yco2 * sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O) * (yh2o - yco2) +
                   xh2o * xco2 * sqrt(aH2O * aCO2) * (kCO2H2O - kH2OCO2)) /
                  amix;
 
   lnPhiCO2 += (term3 - bCO2 / bmix) * amix / (bmix * _Rbar * t15) * log(V / (V + bmix));
 
   return exp(lnPhiCO2);
 }

◆ fugacityCoefficientH2OHighTemp()

ADReal PorousFlowBrineCO2::fugacityCoefficientH2OHighTemp	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	co2_density,
		const ADReal &	xco2,
		const ADReal &	yh2o
	)		const

Definition at line 458 of file PorousFlowBrineCO2.C.

Referenced by fugacityCoefficientsHighTemp().

 {
   using std::sqrt, std::pow, std::log, std::exp;
 
   // Need pressure in bar
   const ADReal pbar = pressure * 1.0e-5;
   // Molar volume in cm^3/mol
   const ADReal V = _Mco2 / co2_density * 1.0e6;
 
   // Redlich-Kwong parameters
   const ADReal yco2 = 1.0 - yh2o;
   const ADReal xh2o = 1.0 - xco2;
 
   const ADReal aCO2 = 8.008e7 - 4.984e4 * temperature;
   const ADReal aH2O = 1.337e8 - 1.4e4 * temperature;
   const Real bCO2 = 28.25;
   const Real bH2O = 15.7;
   const ADReal KH2OCO2 = 1.427e-2 - 4.037e-4 * temperature;
   const ADReal KCO2H2O = 0.4228 - 7.422e-4 * temperature;
   const ADReal kH2OCO2 = KH2OCO2 * yh2o + KCO2H2O * yco2;
   const ADReal kCO2H2O = KCO2H2O * yh2o + KH2OCO2 * yco2;
 
   const ADReal aH2OCO2 = sqrt(aCO2 * aH2O) * (1.0 - kH2OCO2);
   const ADReal aCO2H2O = sqrt(aCO2 * aH2O) * (1.0 - kCO2H2O);
 
   const ADReal amix = yh2o * yh2o * aH2O + yh2o * yco2 * (aH2OCO2 + aCO2H2O) + yco2 * yco2 * aCO2;
   const ADReal bmix = yh2o * bH2O + yco2 * bCO2;
 
   const ADReal t15 = pow(temperature, 1.5);
 
   ADReal lnPhiH2O = bH2O / bmix * (pbar * V / (_Rbar * temperature) - 1.0) -
                     log(pbar * (V - bmix) / (_Rbar * temperature));
   ADReal term3 = (2.0 * yh2o * aH2O + yco2 * (aH2OCO2 + aCO2H2O) -
                   yh2o * yco2 * sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O) * (yh2o - yco2) +
                   xh2o * xco2 * sqrt(aH2O * aCO2) * (kH2OCO2 - kCO2H2O)) /
                  amix;
   term3 -= bH2O / bmix;
   term3 *= amix / (bmix * _Rbar * t15) * log(V / (V + bmix));
   lnPhiH2O += term3;
 
   return exp(lnPhiH2O);
 }

◆ fugacityCoefficientsHighTemp()

void PorousFlowBrineCO2::fugacityCoefficientsHighTemp	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	co2_density,
		const ADReal &	xco2,
		const ADReal &	yh2o,
		ADReal &	fco2,
		ADReal &	fh2o
	)		const

Fugacity coefficients for H2O and CO2 at elevated temperatures (100C < T <= 300C).

Eq. (A8) from Spycher & Pruess (2010)

Parameters

	pressure	gas pressure (Pa)
	temperature	temperature (K)
	co2_density	CO2 density (kg/m^3)
	xco2	mole fraction of CO2 in liquid phase (-)
	yh2o	mole fraction of H2O in gas phase (-)
[out]	fco2	fugacity coefficient for CO2
[out]	fh2o	fugacity coefficient for H2O

Definition at line 440 of file PorousFlowBrineCO2.C.

Referenced by funcABHighTemp().

 {
   if (temperature.value() <= 373.15)
     mooseError(name(),
                ": fugacityCoefficientsHighTemp() is not valid for T <= 373.15K. Use "
                "fugacityCoefficientsLowTemp() instead");
 
   fh2o = fugacityCoefficientH2OHighTemp(pressure, temperature, co2_density, xco2, yh2o);
   fco2 = fugacityCoefficientCO2HighTemp(pressure, temperature, co2_density, xco2, yh2o);
 }

◆ fugacityCoefficientsLowTemp()

void PorousFlowBrineCO2::fugacityCoefficientsLowTemp	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	co2_density,
		ADReal &	fco2,
		ADReal &	fh2o
	)		const

Fugacity coefficients for H2O and CO2 for T <= 100C Eq.

(B7) from Spycher, Pruess and Ennis-King (2003)

Parameters

	pressure	gas pressure (Pa)
	temperature	temperature (K)
	co2_density	CO2 density (kg/m^3)
[out]	fco2	fugacity coefficient for CO2
[out]	fh2o	fugacity coefficient for H2O

Definition at line 397 of file PorousFlowBrineCO2.C.

Referenced by funcABLowTemp().

 {
   using std::log, std::exp, std::pow;
 
   if (temperature.value() > 373.15)
     mooseError(name(),
                ": fugacityCoefficientsLowTemp() is not valid for T > 373.15K. Use "
                "fugacityCoefficientsHighTemp() instead");
 
   // Need pressure in bar
   const ADReal pbar = pressure * 1.0e-5;
 
   // Molar volume in cm^3/mol
   const ADReal V = _Mco2 / co2_density * 1.0e6;
 
   // Redlich-Kwong parameters
   const ADReal aCO2 = 7.54e7 - 4.13e4 * temperature;
   const Real bCO2 = 27.8;
   const Real aCO2H2O = 7.89e7;
   const Real bH2O = 18.18;
 
   const ADReal t15 = pow(temperature, 1.5);
 
   // The fugacity coefficients for H2O and CO2
   auto lnPhi = [V, aCO2, bCO2, t15, this](ADReal a, ADReal b)
   {
     return log(V / (V - bCO2)) + b / (V - bCO2) -
            2.0 * a / (_Rbar * t15 * bCO2) * log((V + bCO2) / V) +
            aCO2 * b / (_Rbar * t15 * bCO2 * bCO2) * (log((V + bCO2) / V) - bCO2 / (V + bCO2));
   };
 
   const ADReal lnPhiH2O = lnPhi(aCO2H2O, bH2O) - log(pbar * V / (_Rbar * temperature));
   const ADReal lnPhiCO2 = lnPhi(aCO2, bCO2) - log(pbar * V / (_Rbar * temperature));
 
   fh2o = exp(lnPhiH2O);
   fco2 = exp(lnPhiCO2);
 }

◆ funcABHighTemp() [1/2]

void PorousFlowBrineCO2::funcABHighTemp	(	Real	pressure,
		Real	temperature,
		Real	Xnacl,
		Real	co2_density,
		Real	xco2,
		Real	yh2o,
		Real &	A,
		Real &	B
	)		const

protected

The function A (Eq.

(11) of Spycher, Pruess and Ennis-King (2003) for T <= 100C, and Eqs. (10) and (17) of Spycher and Pruess (2010) for T > 100C)

Parameters

	pressure	gas pressure (Pa)
	temperature	fluid temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
	co2_density	CO2 density (kg/m^3)
	xco2	mole fraction of CO2 in liquid phase (-)
	yh2o	mole fraction of H2O in gas phase (-)
[out]	A	the function A
[out]	B	the function B

Definition at line 882 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMoleFractions(), funcABHighTemp(), and solveEquilibriumMoleFractionHighTemp().

 {
   using std::exp;
 
   if (temperature <= 373.15)
     mooseError(name(),
                ": funcABHighTemp() is not valid for T <= 373.15K. Use funcABLowTemp() instead");
 
   // Pressure in bar
   const Real pbar = pressure * 1.0e-5;
   // Temperature in C
   const Real Tc = temperature - _T_c2k;
 
   // Reference pressure and partial molar volumes
   const Real pref = -1.9906e-1 + 2.0471e-3 * Tc + 1.0152e-4 * Tc * Tc - 1.4234e-6 * Tc * Tc * Tc +
                     1.4168e-8 * Tc * Tc * Tc * Tc;
   const Real vCO2 = 32.6 + 3.413e-2 * (Tc - 100.0);
   const Real vH2O = 18.1 + 3.137e-2 * (Tc - 100.0);
 
   const Real delta_pbar = pbar - pref;
   const Real Rt = _Rbar * temperature;
 
   // Equilibrium constants
   // Use dummy ADReal temperature as derivatives aren't required
   const ADReal T = temperature;
   Real K0H2O = equilibriumConstantH2O(T).value();
   Real K0CO2 = equilibriumConstantCO2(T).value();
 
   // Fugacity coefficients
   // Use dummy ADReal variables as derivatives aren't required
   const ADReal p = pressure;
   const ADReal rhoco2 = co2_density;
   const ADReal x = xco2;
   const ADReal y = yh2o;
 
   ADReal phiH2O, phiCO2;
   fugacityCoefficientsHighTemp(p, T, rhoco2, x, y, phiCO2, phiH2O);
 
   // Activity coefficients
   const Real gammaH2O = activityCoefficientH2O(T, x).value();
   const Real gammaCO2 = activityCoefficientCO2(T, x).value();
 
   // Activity coefficient for CO2 in brine
   // Use dummy ADReal Xnacl as derivatives aren't required
   const ADReal X = Xnacl;
   const Real gamma = activityCoefficientHighTemp(T, X).value();
 
   A = K0H2O * gammaH2O / (phiH2O.value() * pbar) * exp(delta_pbar * vH2O / Rt);
   B = phiCO2.value() * pbar / (_invMh2o * K0CO2 * gamma * gammaCO2) * exp(-delta_pbar * vCO2 / Rt);
 }

◆ funcABHighTemp() [2/2]

void PorousFlowBrineCO2::funcABHighTemp	(	Real	pressure,
		Real	temperature,
		Real	Xnacl,
		Real	co2_density,
		Real	xco2,
		Real	yh2o,
		Real &	A,
		Real &	dA_dp,
		Real &	dA_dT,
		Real &	B,
		Real &	dB_dp,
		Real &	dB_dT,
		Real &	dB_dX
	)		const

protected

Definition at line 941 of file PorousFlowBrineCO2.C.

 {
   funcABHighTemp(pressure, temperature, Xnacl, co2_density, xco2, yh2o, A, B);
 
   // Use finite differences for derivatives in the high temperature regime
   const Real dp = 1.0e-2;
   const Real dT = 1.0e-6;
   const Real dX = 1.0e-8;
 
   Real A2, B2;
   funcABHighTemp(pressure + dp, temperature, Xnacl, co2_density, xco2, yh2o, A2, B2);
   dA_dp = (A2 - A) / dp;
   dB_dp = (B2 - B) / dp;
 
   funcABHighTemp(pressure, temperature + dT, Xnacl, co2_density, xco2, yh2o, A2, B2);
   dA_dT = (A2 - A) / dT;
   dB_dT = (B2 - B) / dT;
 
   funcABHighTemp(pressure, temperature, Xnacl + dX, co2_density, xco2, yh2o, A2, B2);
   dB_dX = (B2 - B) / dX;
 }

◆ funcABLowTemp()

void PorousFlowBrineCO2::funcABLowTemp	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	co2_density,
		ADReal &	A,
		ADReal &	B
	)		const

protected

Definition at line 846 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMoleFractionsLowTemp().

 {
   using std::exp;
 
   if (temperature.value() > 373.15)
     mooseError(name(),
                ": funcABLowTemp() is not valid for T > 373.15K. Use funcABHighTemp() instead");
 
   // Pressure in bar
   const ADReal pbar = pressure * 1.0e-5;
 
   // Reference pressure and partial molar volumes
   const Real pref = 1.0;
   const Real vCO2 = 32.6;
   const Real vH2O = 18.1;
 
   const ADReal delta_pbar = pbar - pref;
   const ADReal Rt = _Rbar * temperature;
 
   // Equilibrium constants
   const ADReal K0H2O = equilibriumConstantH2O(temperature);
   const ADReal K0CO2 = equilibriumConstantCO2(temperature);
 
   // Fugacity coefficients
   ADReal phiH2O, phiCO2;
   fugacityCoefficientsLowTemp(pressure, temperature, co2_density, phiCO2, phiH2O);
 
   A = K0H2O / (phiH2O * pbar) * exp(delta_pbar * vH2O / Rt);
   B = phiCO2 * pbar / (_invMh2o * K0CO2) * exp(-delta_pbar * vCO2 / Rt);
 }

◆ gasComponentIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::gasComponentIndex ( ) const

inlineinherited

The index of the gas fluid component.

Returns: gas fluid component number

Definition at line 35 of file PorousFlowFluidStateMultiComponentBase.h.

35 { return _gas_fluid_component; };

PorousFlowFluidStateMultiComponentBase::_gas_fluid_component

unsigned int _gas_fluid_component

Fluid component number of the gas phase.

Definition: PorousFlowFluidStateMultiComponentBase.h:101

◆ gasPhaseIndex()

unsigned int PorousFlowFluidStateBase::gasPhaseIndex ( ) const

inlineinherited

The index of the gas phase.

Returns: gas phase number

Definition at line 83 of file PorousFlowFluidStateBase.h.

83 { return _gas_phase_number; };

PorousFlowFluidStateBase::_gas_phase_number

unsigned int _gas_phase_number

Phase number of the gas phase.

Definition: PorousFlowFluidStateBase.h:104

◆ gasProperties()

void PorousFlowBrineCO2::gasProperties	(	const ADReal &	pressure,
		const ADReal &	temperature,
		std::vector< FluidStateProperties > &	fsp
	)		const

Thermophysical properties of the gaseous state.

Parameters

	pressure	gas pressure (Pa)
	temperature	temperature (K)
[out]	FluidStateProperties	data structure

Definition at line 242 of file PorousFlowBrineCO2.C.

Referenced by thermophysicalProperties(), totalMassFraction(), and twoPhaseProperties().

 {
   FluidStateProperties & gas = fsp[_gas_phase_number];
 
   // Gas density, viscosity and enthalpy are approximated with pure CO2 - no correction due
   // to the small amount of water vapor is made
   ADReal co2_density, co2_viscosity;
   _co2_fp.rho_mu_from_p_T(pressure, temperature, co2_density, co2_viscosity);
 
   ADReal co2_enthalpy = _co2_fp.h_from_p_T(pressure, temperature);
 
   // Save the values to the FluidStateProperties object. Note that derivatives wrt z are 0
   gas.density = co2_density;
   gas.viscosity = co2_viscosity;
   gas.enthalpy = co2_enthalpy;
 
   mooseAssert(gas.density.value() > 0.0, "Gas density must be greater than zero");
   gas.internal_energy = gas.enthalpy - pressure / gas.density;
 }

◆ getPressureIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getPressureIndex ( ) const

inlineinherited

Definition at line 92 of file PorousFlowFluidStateMultiComponentBase.h.

Referenced by PorousFlowWaterNCGTest::buildObjects(), and PorousFlowBrineCO2Test::buildObjects().

92 { return _pidx; };

PorousFlowFluidStateMultiComponentBase::_pidx

const unsigned int _pidx

Index of derivative wrt pressure.

Definition: PorousFlowFluidStateMultiComponentBase.h:105

◆ getTemperatureIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getTemperatureIndex ( ) const

inlineinherited

Definition at line 93 of file PorousFlowFluidStateMultiComponentBase.h.

Referenced by PorousFlowWaterNCGTest::buildObjects(), and PorousFlowBrineCO2Test::buildObjects().

93 { return _Tidx; };

PorousFlowFluidStateMultiComponentBase::_Tidx

const unsigned int _Tidx

Index of derivative wrt temperature.

Definition: PorousFlowFluidStateMultiComponentBase.h:109

◆ getXIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getXIndex ( ) const

inlineinherited

Definition at line 95 of file PorousFlowFluidStateMultiComponentBase.h.

Referenced by PorousFlowBrineCO2Test::buildObjects().

95 { return _Xidx; };

PorousFlowFluidStateMultiComponentBase::_Xidx

const unsigned int _Xidx

Index of derivative wrt salt mass fraction X.

Definition: PorousFlowFluidStateMultiComponentBase.h:111

◆ getZIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getZIndex ( ) const

inlineinherited

Definition at line 94 of file PorousFlowFluidStateMultiComponentBase.h.

Referenced by PorousFlowWaterNCGTest::buildObjects(), and PorousFlowBrineCO2Test::buildObjects().

94 { return _Zidx; };

PorousFlowFluidStateMultiComponentBase::_Zidx

const unsigned int _Zidx

Index of derivative wrt total mass fraction Z.

Definition: PorousFlowFluidStateMultiComponentBase.h:107

◆ henryConstant()

ADReal PorousFlowBrineCO2::henryConstant	(	const ADReal &	temperature,
		const ADReal &	Xnacl
	)		const

Henry's constant of dissolution of gas phase CO2 in brine.

From Battistelli et al, A fluid property module for the TOUGH2 simulator for saline brines with non-condensible gas, Proc. Eighteenth Workshop on Geothermal Reservoir Engineering (1993)

Parameters

temperature	fluid temperature (K)
Xnacl	NaCl mass fraction (kg/kg)

Returns: Henry's constant (Pa)

Definition at line 1097 of file PorousFlowBrineCO2.C.

Referenced by enthalpyOfDissolutionGas().

 {
   using std::pow;
 
   // Henry's constant for dissolution in water
   const ADReal Kh_h2o = _brine_fp.henryConstant(temperature, _co2_henry);
 
   // The correction to salt is obtained through the salting out coefficient
   const std::vector<Real> b{1.19784e-1, -7.17823e-4, 4.93854e-6, -1.03826e-8, 1.08233e-11};
 
   // Need temperature in Celsius
   const ADReal Tc = temperature - _T_c2k;
 
   ADReal kb = 0.0;
   for (unsigned int i = 0; i < b.size(); ++i)
     kb += b[i] * pow(Tc, i);
 
   // Need salt mass fraction in molality
   const ADReal xmol = Xnacl / (1.0 - Xnacl) / _Mnacl;
 
   // Henry's constant and its derivative wrt temperature and salt mass fraction
   return Kh_h2o * pow(10.0, xmol * kb);
 }

◆ initialize()

void PorousFlowFluidStateBase::initialize ( )

inlinefinalvirtualinherited

Implements GeneralUserObject.

Definition at line 57 of file PorousFlowFluidStateBase.h.

57 {};

◆ liquidProperties()

void PorousFlowBrineCO2::liquidProperties	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		std::vector< FluidStateProperties > &	fsp
	)		const

Thermophysical properties of the liquid state.

Parameters

	pressure	liquid pressure (Pa)
	temperature	temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
[out]	FluidStateProperties	data structure

Definition at line 265 of file PorousFlowBrineCO2.C.

Referenced by thermophysicalProperties(), totalMassFraction(), and twoPhaseProperties().

 {
   FluidStateProperties & liquid = fsp[_aqueous_phase_number];
 
   // The liquid density includes the density increase due to dissolved CO2
   const ADReal brine_density = _brine_fp.rho_from_p_T_X(pressure, temperature, Xnacl);
 
   // Mass fraction of CO2 in liquid phase
   const ADReal Xco2 = liquid.mass_fraction[_gas_fluid_component];
 
   // The liquid density
   const ADReal co2_partial_density = partialDensityCO2(temperature);
 
   const ADReal liquid_density = 1.0 / (Xco2 / co2_partial_density + (1.0 - Xco2) / brine_density);
 
   // Assume that liquid viscosity is just the brine viscosity
   const ADReal liquid_viscosity = _brine_fp.mu_from_p_T_X(pressure, temperature, Xnacl);
 
   // Liquid enthalpy (including contribution due to the enthalpy of dissolution)
   const ADReal brine_enthalpy = _brine_fp.h_from_p_T_X(pressure, temperature, Xnacl);
 
   // Enthalpy of CO2
   const ADReal co2_enthalpy = _co2_fp.h_from_p_T(pressure, temperature);
 
   // Enthalpy of dissolution
   const ADReal hdis = enthalpyOfDissolution(temperature);
 
   const ADReal liquid_enthalpy = (1.0 - Xco2) * brine_enthalpy + Xco2 * (co2_enthalpy + hdis);
 
   // Save the values to the FluidStateProperties object
   liquid.density = liquid_density;
   liquid.viscosity = liquid_viscosity;
   liquid.enthalpy = liquid_enthalpy;
 
   mooseAssert(liquid.density.value() > 0.0, "Liquid density must be greater than zero");
   liquid.internal_energy = liquid.enthalpy - pressure / liquid.density;
 }

◆ massFractions()

void PorousFlowBrineCO2::massFractions	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		const ADReal &	Z,
		FluidStatePhaseEnum &	phase_state,
		std::vector< FluidStateProperties > &	fsp
	)		const

Mass fractions of CO2 and H2O in both phases, as well as derivatives wrt PorousFlow variables.

Values depend on the phase state (liquid, gas or two phase)

Parameters

	pressure	phase pressure (Pa)
	temperature	phase temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
	Z	total mass fraction of CO2 component
[out]	PhaseStateEnum	current phase state
[out]	FluidStateProperties	data structure

Definition at line 171 of file PorousFlowBrineCO2.C.

Referenced by thermophysicalProperties().

 {
   FluidStateProperties & liquid = fsp[_aqueous_phase_number];
   FluidStateProperties & gas = fsp[_gas_phase_number];
 
   ADReal Xco2 = 0.0;
   ADReal Yh2o = 0.0;
   ADReal Yco2 = 0.0;
 
   // If the amount of CO2 is less than the smallest solubility, then all CO2 will
   // be dissolved, and the equilibrium mass fractions do not need to be computed
   if (Z.value() < _Zmin)
     phase_state = FluidStatePhaseEnum::LIQUID;
 
   else
   {
     // Equilibrium mass fraction of CO2 in liquid and H2O in gas phases
     equilibriumMassFractions(pressure, temperature, Xnacl, Xco2, Yh2o);
 
     Yco2 = 1.0 - Yh2o;
 
     // Determine which phases are present based on the value of z
     phaseState(Z.value(), Xco2.value(), Yco2.value(), phase_state);
   }
 
   // The equilibrium mass fractions calculated above are only correct in the two phase
   // state. If only liquid or gas phases are present, the mass fractions are given by
   // the total mass fraction z
   ADReal Xh2o = 0.0;
 
   switch (phase_state)
   {
     case FluidStatePhaseEnum::LIQUID:
     {
       Xco2 = Z;
       Yco2 = 0.0;
       Xh2o = 1.0 - Z;
       Yh2o = 0.0;
       break;
     }
 
     case FluidStatePhaseEnum::GAS:
     {
       Xco2 = 0.0;
       Yco2 = Z;
       Yh2o = 1.0 - Z;
       break;
     }
 
     case FluidStatePhaseEnum::TWOPHASE:
     {
       // Keep equilibrium mass fractions
       Xh2o = 1.0 - Xco2;
       break;
     }
   }
 
   // Save the mass fractions in the FluidStateProperties object
   liquid.mass_fraction[_aqueous_fluid_component] = Xh2o;
   liquid.mass_fraction[_gas_fluid_component] = Xco2;
   liquid.mass_fraction[_salt_component] = Xnacl;
   gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
   gas.mass_fraction[_gas_fluid_component] = Yco2;
 }

◆ numComponents()

unsigned int PorousFlowFluidStateBase::numComponents ( ) const

inlineinherited

The maximum number of components in this model.

Returns: number of components

Definition at line 71 of file PorousFlowFluidStateBase.h.

71 { return _num_components; };

PorousFlowFluidStateBase::_num_components

unsigned int _num_components

Number of components.

Definition: PorousFlowFluidStateBase.h:100

◆ numPhases()

unsigned int PorousFlowFluidStateBase::numPhases ( ) const

inlineinherited

The maximum number of phases in this model.

Returns: number of phases

Definition at line 65 of file PorousFlowFluidStateBase.h.

Referenced by PorousFlowFluidStateSingleComponentTempl< is_ad >::PorousFlowFluidStateSingleComponentTempl(), and PorousFlowFluidStateTempl< is_ad >::PorousFlowFluidStateTempl().

65 { return _num_phases; };

PorousFlowFluidStateBase::_num_phases

unsigned int _num_phases

Number of phases.

Definition: PorousFlowFluidStateBase.h:98

◆ partialDensityCO2()

ADReal PorousFlowBrineCO2::partialDensityCO2 ( const ADReal & temperature ) const

Partial density of dissolved CO2 From Garcia, Density of aqueous solutions of CO2, LBNL-49023 (2001)

Parameters

temperature fluid temperature (K)

Returns: partial molar density (kg/m^3)

Definition at line 1043 of file PorousFlowBrineCO2.C.

Referenced by liquidProperties(), and saturation().

 {
   // This correlation uses temperature in C
   const ADReal Tc = temperature - _T_c2k;
   // The parial molar volume
   const ADReal V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
 
   return 1.0e6 * _Mco2 / V;
 }

◆ phaseState()

void PorousFlowFluidStateMultiComponentBase::phaseState	(	Real	Zi,
		Real	Xi,
		Real	Yi,
		FluidStatePhaseEnum &	phase_state
	)		const

inherited

Determines the phase state gven the total mass fraction and equilibrium mass fractions.

Parameters

	Zi	total mass fraction
	Xi	equilibrium mass fraction in liquid
	Yi	equilibrium mass fraction in gas
[out]	phase_state	the phase state (gas, liquid, two phase)

Definition at line 36 of file PorousFlowFluidStateMultiComponentBase.C.

Referenced by PorousFlowWaterNCG::massFractions(), and massFractions().

 {
   if (Zi <= Xi)
   {
     // In this case, there is not enough component i to form a gas phase,
     // so only a liquid phase is present
     phase_state = FluidStatePhaseEnum::LIQUID;
   }
   else if (Zi > Xi && Zi < Yi)
   {
     // Two phases are present
     phase_state = FluidStatePhaseEnum::TWOPHASE;
   }
   else // (Zi >= Yi)
   {
     // In this case, there is not enough water to form a liquid
     // phase, so only a gas phase is present
     phase_state = FluidStatePhaseEnum::GAS;
   }
 }

◆ rachfordRice()

Real PorousFlowFluidStateFlash::rachfordRice	(	Real	vf,
		std::vector< Real > &	Zi,
		std::vector< Real > &	Ki
	)		const

inherited

Rachford-Rice equation for vapor fraction.

Can be solved analytically for two components in two phases, but must be solved iteratively using a root finding algorithm for more components. This equation has the nice property that it is monotonic in the interval [0,1], so that only a small number of iterations are typically required to find the root.

The Rachford-Rice equation can also be used to check whether the phase state is two phase, single phase gas, or single phase liquid. Evaluate f(v), the Rachford-Rice equation evaluated at the vapor mass fraction.

If f(0) < 0, then the mixture is below the bubble point, and only a single phase liquid can exist

If f(1) > 0, then the mixture is above the dew point, and only a single phase gas exists.

If f(0) >= 0 and f(1) <= 0, the mixture is between the bubble and dew points, and both gas and liquid phases exist.

Parameters

vf	vapor fraction
Zi	mass fractions
Ki	equilibrium constants

Returns: f(x)

Definition at line 26 of file PorousFlowFluidStateFlash.C.

Referenced by PorousFlowFluidStateFlash::vaporMassFraction().

 {
   const std::size_t num_z = Zi.size();
   // Check that the sizes of the mass fractions and equilibrium constant vectors are correct
   if (Ki.size() != num_z + 1)
     mooseError("The number of mass fractions or equilibrium components passed to rachfordRice is "
                "not correct");
 
   Real f = 0.0;
   Real Z_total = 0.0;
 
   for (std::size_t i = 0; i < num_z; ++i)
   {
     f += Zi[i] * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0));
     Z_total += Zi[i];
   }
 
   // Add the last component (with total mass fraction = 1 - z_total)
   f += (1.0 - Z_total) * (Ki[num_z] - 1.0) / (1.0 + x * (Ki[num_z] - 1.0));
 
   return f;
 }

◆ rachfordRiceDeriv()

Real PorousFlowFluidStateFlash::rachfordRiceDeriv	(	Real	vf,
		std::vector< Real > &	Zi,
		std::vector< Real > &	Ki
	)		const

inherited

Derivative of Rachford-Rice equation wrt vapor fraction.

Has the nice property that it is strictly negative in the interval [0,1]

Parameters

vf	vapor fraction
Zi	mass fractions
Ki	equilibrium constants

Returns: f'(x)

Definition at line 52 of file PorousFlowFluidStateFlash.C.

Referenced by PorousFlowFluidStateFlash::vaporMassFraction().

 {
   const std::size_t num_Z = Zi.size();
   // Check that the sizes of the mass fractions and equilibrium constant vectors are correct
   if (Ki.size() != num_Z + 1)
     mooseError("The number of mass fractions or equilibrium components passed to rachfordRice is "
                "not correct");
 
   Real df = 0.0;
   Real Z_total = 0.0;
 
   for (std::size_t i = 0; i < num_Z; ++i)
   {
     df -= Zi[i] * (Ki[i] - 1.0) * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0)) /
           (1.0 + x * (Ki[i] - 1.0));
     Z_total += Zi[i];
   }
 
   // Add the last component (with total mass fraction = 1 - z_total)
   df -= (1.0 - Z_total) * (Ki[num_Z] - 1.0) * (Ki[num_Z] - 1.0) / (1.0 + x * (Ki[num_Z] - 1.0)) /
         (1.0 + x * (Ki[num_Z] - 1.0));
 
   return df;
 }

◆ saltComponentIndex()

unsigned int PorousFlowBrineCO2::saltComponentIndex ( ) const

inline

The index of the salt component.

Returns: salt component number

Definition at line 321 of file PorousFlowBrineCO2.h.

321 { return _salt_component; };

PorousFlowBrineCO2::_salt_component

const unsigned int _salt_component

Salt component index.

Definition: PorousFlowBrineCO2.h:460

◆ saturation()

ADReal PorousFlowBrineCO2::saturation	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		const ADReal &	Z,
		std::vector< FluidStateProperties > &	fsp
	)		const

Gas saturation in the two-phase region.

Parameters

pressure	gas pressure (Pa)
temperature	phase temperature (K)
Xnacl	NaCl mass fraction (kg/kg)
Z	total mass fraction of CO2 component
FluidStateProperties	data structure

Returns: gas saturation (-)

Definition at line 307 of file PorousFlowBrineCO2.C.

Referenced by totalMassFraction(), and twoPhaseProperties().

 {
   auto & gas = fsp[_gas_phase_number];
   auto & liquid = fsp[_aqueous_fluid_component];
 
   // Approximate liquid density as saturation isn't known yet, by using the gas
   // pressure rather than the liquid pressure. This does result in a small error
   // in the calculated saturation, but this is below the error associated with
   // the correlations. A more accurate saturation could be found iteraviely,
   // at the cost of increased computational expense
 
   // Gas density
   const ADReal gas_density = _co2_fp.rho_from_p_T(pressure, temperature);
 
   // Approximate liquid density as saturation isn't known yet
   const ADReal brine_density = _brine_fp.rho_from_p_T_X(pressure, temperature, Xnacl);
 
   // Mass fraction of CO2 in liquid phase
   const ADReal Xco2 = liquid.mass_fraction[_gas_fluid_component];
 
   // The liquid density
   const ADReal co2_partial_density = partialDensityCO2(temperature);
 
   const ADReal liquid_density = 1.0 / (Xco2 / co2_partial_density + (1.0 - Xco2) / brine_density);
 
   const ADReal Yco2 = gas.mass_fraction[_gas_fluid_component];
 
   // Set mass equilibrium constants used in the calculation of vapor mass fraction
   const ADReal K0 = Yco2 / Xco2;
   const ADReal K1 = (1.0 - Yco2) / (1.0 - Xco2);
   const ADReal vapor_mass_fraction = vaporMassFraction(Z, K0, K1);
 
   // The gas saturation in the two phase case
   const ADReal saturation = vapor_mass_fraction * liquid_density /
                             (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
 
   return saturation;
 }

◆ smoothCubicInterpolation()

void PorousFlowBrineCO2::smoothCubicInterpolation	(	Real	temperature,
		Real	f0,
		Real	df0,
		Real	f1,
		Real	df1,
		Real &	value,
		Real &	deriv
	)		const

protected

Cubic function to smoothly interpolate between the low temperature and elevated temperature models for 99C < T < 109C.

Parameters

	temperature	temperature (K)
	f0	function value at T = 372K (99C)
	df0	derivative of function at T = 372K (99C)
	f1	function value at T = 382K (109C)
	df1	derivative of function at T = 382K (109C)
[out]	value	value at the given temperature
[out]	deriv	derivative at the given temperature

Definition at line 1162 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMoleFractions().

 {
   // Coefficients of cubic polynomial
   const Real dT = _Tupper - _Tlower;
 
   const Real a = f0;
   const Real b = df0 * dT;
   const Real c = 3.0 * (f1 - f0) - (2.0 * df0 + df1) * dT;
   const Real d = 2.0 * (f0 - f1) + (df0 + df1) * dT;
 
   const Real t2 = temperature * temperature;
   const Real t3 = temperature * t2;
 
   value = a + b * temperature + c * t2 + d * t3;
   deriv = b + 2.0 * c * temperature + 3.0 * d * t2;
 }

◆ solveEquilibriumMoleFractionHighTemp()

void PorousFlowBrineCO2::solveEquilibriumMoleFractionHighTemp	(	Real	pressure,
		Real	temperature,
		Real	Xnacl,
		Real	co2_density,
		Real &	xco2,
		Real &	yh2o
	)		const

Function to solve for yh2o and xco2 iteratively in the elevated temperature regime (T > 100C)

Parameters

	pressure	gas pressure (Pa)
	temperature	fluid temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
	co2_density	CO2 density (kg/m^3)
[out]	xco2	mole fraction of CO2 in liquid phase (-)
[out]	yh2o	mole fraction of H2O in gas phase (-)

Definition at line 976 of file PorousFlowBrineCO2.C.

Referenced by equilibriumMoleFractions().

 {
   // Initial guess for yh2o and xco2 (from Spycher and Pruess (2010))
   Real y = _brine_fp.vaporPressure(temperature, 0.0) / pressure;
   Real x = 0.009;
 
   // Need salt mass fraction in molality
   const Real mnacl = Xnacl / (1.0 - Xnacl) / _Mnacl;
 
   // If y > 1, then just use y = 1, x = 0 (only a gas phase)
   if (y >= 1.0)
   {
     y = 1.0;
     x = 0.0;
   }
   else
   {
     // Residual function for Newton-Raphson
     auto fy = [mnacl, this](Real y, Real A, Real B) {
       return y -
              (1.0 - B) * _invMh2o / ((1.0 / A - B) * (2.0 * mnacl + _invMh2o) + 2.0 * mnacl * B);
     };
 
     // Derivative of fy wrt y
     auto dfy = [mnacl, this](Real A, Real B, Real dA, Real dB)
     {
       const Real denominator = (1.0 / A - B) * (2.0 * mnacl + _invMh2o) + 2.0 * mnacl * B;
       return 1.0 + _invMh2o * dB / denominator +
              (1.0 - B) * _invMh2o *
                  (2.0 * mnacl * dB - (2.0 * mnacl + _invMh2o) * (dB + dA / A / A)) / denominator /
                  denominator;
     };
 
     Real A, B;
     Real dA, dB;
     const Real dy = 1.0e-8;
 
     // Solve for yh2o using Newton-Raphson method
     unsigned int iter = 0;
     const Real tol = 1.0e-12;
     const unsigned int max_its = 10;
     funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y, A, B);
 
     while (std::abs(fy(y, A, B)) > tol)
     {
       funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y, A, B);
       // Finite difference derivatives of A and B wrt y
       funcABHighTemp(pressure, temperature, Xnacl, co2_density, x, y + dy, dA, dB);
       dA = (dA - A) / dy;
       dB = (dB - B) / dy;
 
       y = y - fy(y, A, B) / dfy(A, B, dA, dB);
 
       x = B * (1.0 - y);
 
       // Break if not converged and just use the value
       if (iter > max_its)
         break;
     }
   }
 
   yh2o = y;
   xco2 = x;
 }

◆ thermophysicalProperties() [1/2]

void PorousFlowBrineCO2::thermophysicalProperties	(	Real	pressure,
		Real	temperature,
		Real	Xnacl,
		Real	Z,
		unsigned int	qp,
		std::vector< FluidStateProperties > &	fsp
	)		const

overridevirtual

Determines the complete thermophysical state of the system for a given set of primary variables.

Parameters

	pressure	gas phase pressure (Pa)
	temperature	fluid temperature (K)
	Xnacl	mass fraction of NaCl
	Z	total mass fraction of fluid component
	qp	quadpoint index
[out]	fsp	the FluidStateProperties struct containing all properties

Implements PorousFlowFluidStateMultiComponentBase.

Definition at line 91 of file PorousFlowBrineCO2.C.

 {
   // Make AD versions of primary variables then call AD thermophysicalProperties()
   ADReal p = pressure;
   Moose::derivInsert(p.derivatives(), _pidx, 1.0);
   ADReal T = temperature;
   Moose::derivInsert(T.derivatives(), _Tidx, 1.0);
   ADReal Zco2 = Z;
   Moose::derivInsert(Zco2.derivatives(), _Zidx, 1.0);
   ADReal X = Xnacl;
   Moose::derivInsert(X.derivatives(), _Xidx, 1.0);
 
   thermophysicalProperties(p, T, X, Zco2, qp, fsp);
 }

◆ thermophysicalProperties() [2/2]

void PorousFlowBrineCO2::thermophysicalProperties	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		const ADReal &	Z,
		unsigned int	qp,
		std::vector< FluidStateProperties > &	fsp
	)		const

overridevirtual

Implements PorousFlowFluidStateMultiComponentBase.

Definition at line 112 of file PorousFlowBrineCO2.C.

 {
   FluidStateProperties & liquid = fsp[_aqueous_phase_number];
   FluidStateProperties & gas = fsp[_gas_phase_number];
 
   // Check whether the input temperature is within the region of validity
   checkVariables(pressure.value(), temperature.value());
 
   // Clear all of the FluidStateProperties data
   clearFluidStateProperties(fsp);
 
   FluidStatePhaseEnum phase_state;
   massFractions(pressure, temperature, Xnacl, Z, phase_state, fsp);
 
   switch (phase_state)
   {
     case FluidStatePhaseEnum::GAS:
     {
       // Set the gas saturations
       gas.saturation = 1.0;
 
       // Calculate gas properties
       gasProperties(pressure, temperature, fsp);
 
       break;
     }
 
     case FluidStatePhaseEnum::LIQUID:
     {
       // Calculate the liquid properties
       const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0, qp);
       liquidProperties(liquid_pressure, temperature, Xnacl, fsp);
 
       break;
     }
 
     case FluidStatePhaseEnum::TWOPHASE:
     {
       // Calculate the gas and liquid properties in the two phase region
       twoPhaseProperties(pressure, temperature, Xnacl, Z, qp, fsp);
 
       break;
     }
   }
 
   // Liquid saturations can now be set
   liquid.saturation = 1.0 - gas.saturation;
 
   // Save pressures to FluidStateProperties object
   gas.pressure = pressure;
   liquid.pressure = pressure - _pc.capillaryPressure(liquid.saturation, qp);
 }

◆ totalMassFraction()

Real PorousFlowBrineCO2::totalMassFraction	(	Real	pressure,
		Real	temperature,
		Real	Xnacl,
		Real	saturation,
		unsigned int	qp
	)		const

overridevirtual

Total mass fraction of fluid component summed over all phases in the two-phase state for a specified gas saturation.

Parameters

pressure	gas pressure (Pa)
temperature	temperature (K)
Xnacl	NaCl mass fraction (kg/kg)
saturation	gas saturation (-)
qp	quadpoint index

Returns: total mass fraction Z (-)

Implements PorousFlowFluidStateMultiComponentBase.

Definition at line 1054 of file PorousFlowBrineCO2.C.

 {
   // Check whether the input pressure and temperature are within the region of validity
   checkVariables(pressure, temperature);
 
   // As we do not require derivatives, we can simply ignore their initialisation
   const ADReal p = pressure;
   const ADReal T = temperature;
   const ADReal X = Xnacl;
 
   // FluidStateProperties data structure
   std::vector<FluidStateProperties> fsp(_num_phases, FluidStateProperties(_num_components));
   FluidStateProperties & liquid = fsp[_aqueous_phase_number];
   FluidStateProperties & gas = fsp[_gas_phase_number];
 
   // Calculate equilibrium mass fractions in the two-phase state
   ADReal Xco2, Yh2o;
   equilibriumMassFractions(p, T, X, Xco2, Yh2o);
 
   // Save the mass fractions in the FluidStateMassFractions object
   const ADReal Yco2 = 1.0 - Yh2o;
   liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - Xco2;
   liquid.mass_fraction[_gas_fluid_component] = Xco2;
   gas.mass_fraction[_aqueous_fluid_component] = Yh2o;
   gas.mass_fraction[_gas_fluid_component] = Yco2;
 
   // Gas properties
   gasProperties(pressure, temperature, fsp);
 
   // Liquid properties
   const ADReal liquid_saturation = 1.0 - saturation;
   const ADReal liquid_pressure = p - _pc.capillaryPressure(liquid_saturation, qp);
   liquidProperties(liquid_pressure, T, X, fsp);
 
   // The total mass fraction of ncg (z) can now be calculated
   const ADReal Z = (saturation * gas.density * Yco2 + liquid_saturation * liquid.density * Xco2) /
                    (saturation * gas.density + liquid_saturation * liquid.density);
 
   return Z.value();
 }

◆ twoPhaseProperties()

void PorousFlowBrineCO2::twoPhaseProperties	(	const ADReal &	pressure,
		const ADReal &	temperature,
		const ADReal &	Xnacl,
		const ADReal &	Z,
		unsigned int	qp,
		std::vector< FluidStateProperties > &	fsp
	)		const

Gas and liquid properties in the two-phase region.

Parameters

	pressure	gas pressure (Pa)
	temperature	phase temperature (K)
	Xnacl	NaCl mass fraction (kg/kg)
	Z	total mass fraction of NCG component
	qp	quadpoint for capillary presssure
[out]	FluidStateProperties	data structure

Definition at line 351 of file PorousFlowBrineCO2.C.

Referenced by thermophysicalProperties().

 {
   auto & gas = fsp[_gas_phase_number];
 
   // Calculate all of the gas phase properties, as these don't depend on saturation
   gasProperties(pressure, temperature, fsp);
 
   // The gas saturation in the two phase case
   gas.saturation = saturation(pressure, temperature, Xnacl, Z, fsp);
 
   // The liquid pressure and properties can now be calculated
   const ADReal liquid_pressure = pressure - _pc.capillaryPressure(1.0 - gas.saturation, qp);
   liquidProperties(liquid_pressure, temperature, Xnacl, fsp);
 }

◆ validParams()

InputParameters PorousFlowBrineCO2::validParams ( )

static

Definition at line 19 of file PorousFlowBrineCO2.C.

 {
   InputParameters params = PorousFlowFluidStateMultiComponentBase::validParams();
   params.addRequiredParam<UserObjectName>("brine_fp", "The name of the user object for brine");
   params.addRequiredParam<UserObjectName>("co2_fp", "The name of the user object for CO2");
   params.addParam<unsigned int>("salt_component", 2, "The component number of salt");
   params.addClassDescription("Fluid state class for brine and CO2");
   return params;
 }

◆ vaporMassFraction() [1/3]

Real PorousFlowFluidStateFlash::vaporMassFraction	(	Real	Z0,
		Real	K0,
		Real	K1
	)		const

inherited

Solves Rachford-Rice equation to provide vapor mass fraction.

For two components, the analytical solution is used, while for cases with more than two components, a Newton-Raphson iterative solution is calculated.

Parameters

Zi	total mass fraction(s)
Ki	equilibrium constant(s)

Returns: vapor mass fraction

Definition at line 80 of file PorousFlowFluidStateFlash.C.

Referenced by PorousFlowWaterNCG::saturation(), saturation(), and PorousFlowFluidStateFlash::vaporMassFraction().

 {
   return (Z0 * (K1 - K0) - (K1 - 1.0)) / ((K0 - 1.0) * (K1 - 1.0));
 }

◆ vaporMassFraction() [2/3]

ADReal PorousFlowFluidStateFlash::vaporMassFraction	(	const ADReal &	Z0,
		const ADReal &	K0,
		const ADReal &	K1
	)		const

inherited

Definition at line 86 of file PorousFlowFluidStateFlash.C.

 {
   return (Z0 * (K1 - K0) - (K1 - 1.0)) / ((K0 - 1.0) * (K1 - 1.0));
 }

◆ vaporMassFraction() [3/3]

Real PorousFlowFluidStateFlash::vaporMassFraction	(	std::vector< Real > &	Zi,
		std::vector< Real > &	Ki
	)		const

inherited

Definition at line 94 of file PorousFlowFluidStateFlash.C.

 {
   // Check that the sizes of the mass fractions and equilibrium constant vectors are correct
   if (Ki.size() != Zi.size() + 1)
     mooseError("The number of mass fractions or equilibrium components passed to rachfordRice is "
                "not correct");
   Real v;
 
   // If there are only two components, an analytical solution is possible
   if (Ki.size() == 2)
     v = vaporMassFraction(Zi[0], Ki[0], Ki[1]);
   else
   {
     // More than two components - solve the Rachford-Rice equation using
     // Newton-Raphson method
     // Initial guess for vapor mass fraction
     Real v0 = 0.5;
     unsigned int iter = 0;
 
     while (std::abs(rachfordRice(v0, Zi, Ki)) > _nr_tol)
     {
       v0 = v0 - rachfordRice(v0, Zi, Ki) / rachfordRiceDeriv(v0, Zi, Ki);
       iter++;
 
       if (iter > _nr_max_its)
         break;
     }
     v = v0;
   }
   return v;
 }