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.

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.

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.

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.

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.

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

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

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.

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.

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.

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

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

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

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

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.

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.

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

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.

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

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

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.

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).

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).

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)

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

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

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.

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.

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

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

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.

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

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

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

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

virtual Real	spatialValue (const Point &) const

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

void	setPrimaryThreadCopy (UserObject *primary)

UserObject *	primaryThreadCopy ()

SubProblem &	getSubProblem () const

bool	shouldDuplicateInitialExecution () 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)

std::set< UserObjectName >	getDependObjects () const

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

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

unsigned int	systemNumber () const

virtual bool	needThreadedCopy () const

virtual bool	enabled () const

std::shared_ptr< MooseObject >	getSharedPtr ()

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

bool	isKokkosObject () 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	mooseDeprecatedNoTrace (Args &&... args) const

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

void	callMooseError (std::string msg, const bool with_prefix, const hit::Node *node=nullptr, const bool show_trace=true) 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 UserObjectBase &	getUserObjectBase (const std::string &param_name, bool is_dependency=true) const

const UserObjectBase &	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)

std::pair< Moose::Kokkos::MaterialProperty< T, dimension >, std::set< SubdomainID > >	getKokkosBlockMaterialProperty (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

Moose::Kokkos::Function	getKokkosFunction (const std::string &name) const

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

Moose::Kokkos::Function	getKokkosFunctionByName (const FunctionName &name) const

const T &	getKokkosFunctionByName (const FunctionName &name) const

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

bool	hasKokkosFunctionByName (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).

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, const bool show_trace=true)

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, NameFunc &&name_func)

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.

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.

virtual void	addPostprocessorDependencyHelper (const PostprocessorName &name) const override

virtual void	addVectorPostprocessorDependencyHelper (const VectorPostprocessorName &name) const override

virtual void	addUserObjectDependencyHelper (const UserObjectBase &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.

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.

const BrineFluidProperties &	_brine_fp
	Fluid properties UserObject for water.

const SinglePhaseFluidProperties &	_co2_fp
	Fluid properties UserObject for the CO2.

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

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

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

const Real	_Mnacl
	Molar mass of NaCL.

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

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

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

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

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

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

unsigned int	_gas_fluid_component
	Fluid component number of the gas phase.

const unsigned int	_pidx
	Index of derivative wrt pressure.

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

const unsigned int	_Tidx
	Index of derivative wrt temperature.

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

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

const Real	_nr_tol
	Tolerance for Newton-Raphson iterations.

unsigned int	_num_phases
	Number of phases.

unsigned int	_num_components
	Number of components.

const unsigned int	_aqueous_phase_number
	Phase number of the aqueous phase.

unsigned int	_gas_phase_number
	Phase number of the gas phase.

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

const Real	_T_c2k
	Conversion from C to K.

const PorousFlowCapillaryPressure &	_pc
	Capillary pressure UserObject.

FluidStateProperties	_empty_fsp
	Empty FluidStateProperties object.

const Moose::CoordinateSystemType &	_coord_sys

const THREAD_ID	_tid

SubProblem &	_subproblem

FEProblemBase &	_fe_problem

SystemBase &	_sys

Assembly &	_assembly

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

Private Member Functions
const Moose::FunctionBase *	getKokkosFunctionByNameHelper (const FunctionName &name) const

const UserObjectBase &	getUserObjectFromFEProblem (const UserObjectName &object_name, const THREAD_ID tid=0) const

const T &	castUserObject (const UserObjectBase &uo_base, const std::string &param_name="") const

void	mooseObjectError (const std::string &param_name, std::stringstream &oss) const

const std::string &	userObjectType (const UserObjectBase &uo) const

const std::string &	userObjectName (const UserObjectBase &uo) const

const PostprocessorName &	getPostprocessorNameInternal (const std::string &param_name, const unsigned int index, const bool allow_default_value=true) const

bool	isDefaultPostprocessorValueByName (const PostprocessorName &name) const

PostprocessorValue	getDefaultPostprocessorValueByName (const PostprocessorName &name) const

void	checkParam (const std::string &param_name, const unsigned int index=std::numeric_limits< unsigned int >::max()) const

bool	postprocessorsAdded () const

const VectorPostprocessorValue &	getVectorPostprocessorByNameHelper (const VectorPostprocessorName &name, const std::string &vector_name, bool broadcast, std::size_t t_index) const

const VectorPostprocessorContext< VectorPostprocessorValue > &	getVectorPostprocessorContextByNameHelper (const VectorPostprocessorName &name, const std::string &vector_name) const

bool	vectorPostprocessorsAdded () const

bool	reportersAdded () const

void	possiblyCheckHasReporter (const ReporterName &reporter_name, const std::string &param_name="") const

RestartableDataValue &	registerRestartableDataOnApp (std::unique_ptr< RestartableDataValue > data, THREAD_ID tid) const

void	registerRestartableNameWithFilterOnApp (const std::string &name, Moose::RESTARTABLE_FILTER filter)

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

virtual std::string	meshPropertyPrefix (const std::string &data_name) const

const RestartableDataValue &	getMeshPropertyInternal (const std::string &data_name, const std::string &prefix) const

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

const VariableValue *	getDefaultValue (const std::string &var_name) const

const ADVariableValue *	getADDefaultValue (const std::string &var_name) const

void	checkVar (const std::string &var_name) const

void	validateExecutionerType (const std::string &name, const std::string &fn_name) const

Moose::MaterialDataType	getMaterialDataType (const std::set< BoundaryID > &boundary_ids) const

unsigned int	getMaxQps () const

void	addConsumedPropertyName (const MooseObjectName &obj_name, const std::string &prop_name)

const PostprocessorValue &	getPostprocessorValueInternal (const std::string &param_name, unsigned int index, std::size_t t_index) const

const PostprocessorValue &	getPostprocessorValueInternal (const std::string &param_name, unsigned int index, std::size_t t_index) const

const PostprocessorValue &	getPostprocessorValueByNameInternal (const PostprocessorName &name, std::size_t t_index) const

const PostprocessorValue &	getPostprocessorValueByNameInternal (const PostprocessorName &name, std::size_t t_index) const

void	possiblyCheckHasVectorPostprocessor (const std::string &param_name, const std::string &vector_name) const

void	possiblyCheckHasVectorPostprocessor (const std::string &param_name, const std::string &vector_name) const

void	possiblyCheckHasVectorPostprocessorByName (const VectorPostprocessorName &name, const std::string &vector_name) const

void	possiblyCheckHasVectorPostprocessorByName (const VectorPostprocessorName &name, const std::string &vector_name) const

Static Private Member Functions
static const hit::Node *	getHitNode (const InputParameters &params)

static std::string	messagePrefix (const InputParameters &params, const bool hit_prefix)

Private Attributes
UserObject *	_primary_thread_copy

std::set< std::string >	_supplied_uo

const ParallelParamObject &	_parent

const MooseBase &	_si_moose_base

const FEProblemBase *	_si_problem

ExecFlagEnum	_empty_execute_enum

const MooseObject &	_fni_object

const InputParameters &	_fni_params

FEProblemBase &	_fni_feproblem

const THREAD_ID	_fni_tid

const MooseObject &	_uoi_moose_object

const FEProblemBase &	_uoi_feproblem

const THREAD_ID	_uoi_tid

const MooseObject &	_ppi_moose_object

const InputParameters &	_ppi_params

const FEProblemBase &	_ppi_feproblem

std::map< PostprocessorName, std::unique_ptr< PostprocessorValue > >	_default_values

const bool	_broadcast_by_default

const MooseObject &	_vpi_moose_object

const FEProblemBase &	_vpi_feproblem

const THREAD_ID	_vpi_tid

const InputParameters &	_ri_params

FEProblemBase &	_ri_fe_problem_base

const ReporterData &	_ri_reporter_data

const MooseObject &	_ri_moose_object

const InputParameters &	_dni_params

FEProblemBase &	_dni_feproblem

const MooseObject *const	_dni_moose_object_ptr

const InputParameters &	_si_params

FEProblemBase &	_si_feproblem

THREAD_ID	_si_tid

const RestartableDataMapName	_metaname

std::string	_restartable_name

MooseApp &	_meta_data_app

const MooseObject *const	_meta_data_object

const InputParameters &	_sc_parameters

const std::string &	_sc_name

const bool	_sc_is_implicit

std::unordered_map< std::string, std::vector< MooseVariableScalar * > >	_coupled_scalar_vars

std::unordered_map< std::string, std::unique_ptr< VariableValue > >	_default_value

std::unordered_map< std::string, std::unique_ptr< ADVariableValue > >	_dual_default_value

std::vector< MooseVariableScalar * >	_coupled_moose_scalar_vars

std::unordered_map< std::string, std::vector< MooseVariableFieldBase * > >	_sc_coupled_vars

std::set< TagID >	_sc_coupleable_vector_tags

std::set< TagID >	_sc_coupleable_matrix_tags

const MooseObject &	_mi_moose_object

const bool	_mi_boundary_restricted

const std::set< SubdomainID > &	_mi_block_ids

const std::set< BoundaryID > &	_mi_boundary_ids

std::vector< std::unique_ptr< OptionalMaterialPropertyProxyBase< MaterialPropertyInterface > > >	_optional_property_proxies

const std::string	_ti_name

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.

{
  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);
}

Referenced by equilibriumMoleFractionsLowTemp().

◆ 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.

{
  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);
  }
}

Referenced by funcABHighTemp().

◆ 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.

{
  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);
  }
}

Referenced by funcABHighTemp().

◆ 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.

{
  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);
}

Referenced by funcABHighTemp().

◆ 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; };

◆ aqueousPhaseIndex()

unsigned int PorousFlowFluidStateBase::aqueousPhaseIndex ( ) const

inlineinherited

The index of the aqueous phase.

Returns: aqueous phase number

Definition at line 79 of file PorousFlowFluidStateBase.h.

79{ return _aqueous_phase_number; };

◆ checkVariables()

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

protectedvirtual

Check the input variables.

Definition at line 1182 of file PorousFlowBrineCO2.C.

{
  // 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");
}

Referenced by thermophysicalProperties(), and totalMassFraction().

◆ 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.

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

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

◆ 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 1153 of file PorousFlowBrineCO2.C.

{
  // 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;
}

Referenced by liquidProperties().

◆ 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 1123 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.

{
  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);
}

Referenced by funcABHighTemp(), and funcABLowTemp().

◆ 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.

{
  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);
}

Referenced by funcABHighTemp(), and funcABLowTemp().

◆ 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.

{
  // 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;
}

Referenced by massFractions(), and totalMassFraction().

◆ 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.

{
  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);
  }
}

Referenced by equilibriumMassFractions().

◆ 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.

{
  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);
}

Referenced by equilibriumMoleFractions().

◆ execute()

void PorousFlowFluidStateBase::execute ( )

inlinefinalvirtualinherited

Implements GeneralUserObject.

Definition at line 60 of file PorousFlowFluidStateBase.h.

60{};

◆ finalize()

void PorousFlowFluidStateBase::finalize ( )

inlinefinalvirtualinherited

Implements GeneralUserObject.

Definition at line 61 of file PorousFlowFluidStateBase.h.

61{};

◆ 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.

{
  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);
}

Referenced by fugacityCoefficientsHighTemp().

◆ 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.

{
  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);
}

Referenced by fugacityCoefficientsHighTemp().

◆ 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.

{
  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);
}

Referenced by funcABHighTemp().

◆ 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.

{
  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);
}

Referenced by funcABLowTemp().

◆ 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.

{
  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);
}

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

◆ 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.

{
  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);
}

Referenced by equilibriumMoleFractionsLowTemp().

◆ 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; };

◆ gasPhaseIndex()

unsigned int PorousFlowFluidStateBase::gasPhaseIndex ( ) const

inlineinherited

The index of the gas phase.

Returns: gas phase number

Definition at line 85 of file PorousFlowFluidStateBase.h.

85{ return _gas_phase_number; };

◆ 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.

{
  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;
}

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

◆ getPressureIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getPressureIndex ( ) const

inlineinherited

Definition at line 92 of file PorousFlowFluidStateMultiComponentBase.h.

92{ return _pidx; };

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

◆ getTemperatureIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getTemperatureIndex ( ) const

inlineinherited

Definition at line 93 of file PorousFlowFluidStateMultiComponentBase.h.

93{ return _Tidx; };

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

◆ getXIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getXIndex ( ) const

inlineinherited

Definition at line 95 of file PorousFlowFluidStateMultiComponentBase.h.

95{ return _Xidx; };

Referenced by PorousFlowBrineCO2Test::buildObjects().

◆ getZIndex()

unsigned int PorousFlowFluidStateMultiComponentBase::getZIndex ( ) const

inlineinherited

Definition at line 94 of file PorousFlowFluidStateMultiComponentBase.h.

94{ return _Zidx; };

PorousFlowFluidStateMultiComponentBase::_Zidx

const unsigned int _Zidx

Index of derivative wrt total mass fraction Z.

Definition PorousFlowFluidStateMultiComponentBase.h:107

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

◆ 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 1098 of file PorousFlowBrineCO2.C.

{
  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);
}

Referenced by enthalpyOfDissolutionGas().

◆ initialize()

void PorousFlowFluidStateBase::initialize ( )

inlinefinalvirtualinherited

Implements GeneralUserObject.

Definition at line 59 of file PorousFlowFluidStateBase.h.

59{};

◆ 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.

{
  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;
}

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

◆ 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.

{
  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;
}

Referenced by thermophysicalProperties().

◆ numComponents()

unsigned int PorousFlowFluidStateBase::numComponents ( ) const

inlineinherited

The maximum number of components in this model.

Returns: number of components

Definition at line 73 of file PorousFlowFluidStateBase.h.

73{ return _num_components; };

◆ numPhases()

unsigned int PorousFlowFluidStateBase::numPhases ( ) const

inlineinherited

The maximum number of phases in this model.

Returns: number of phases

Definition at line 67 of file PorousFlowFluidStateBase.h.

67{ return _num_phases; };

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

◆ 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 1044 of file PorousFlowBrineCO2.C.

{
  // 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;
}

Referenced by liquidProperties(), and saturation().

◆ 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.

{
  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;
  }
}

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

◆ 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.

{
  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;
}

Referenced by PorousFlowFluidStateFlash::vaporMassFraction().

◆ 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.

{
  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;
}

Referenced by PorousFlowFluidStateFlash::vaporMassFraction().

◆ 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; };

◆ 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.

{
  auto & gas = fsp[_gas_phase_number];
  auto & liquid = fsp[_aqueous_phase_number];
 
  // 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;
}

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

◆ 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 1163 of file PorousFlowBrineCO2.C.

{
  // 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;
}

Referenced by equilibriumMoleFractions().

◆ 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.

{
  // 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;
}

Referenced by equilibriumMoleFractions().

◆ thermophysicalProperties() [1/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);
}

◆ thermophysicalProperties() [2/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);
}

Referenced by thermophysicalProperties().

◆ 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 1055 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.

{
  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);
}

Referenced by thermophysicalProperties().

◆ 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]

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() [2/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.

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

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

◆ 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;
}