Convection of primary species in given equilibrium species. More...

#include <CoupledConvectionReactionSub.h>

Inheritance diagram for CoupledConvectionReactionSub:

Public Types
typedef DerivativeMaterialPropertyNameInterface::SymbolName	SymbolName

Public Member Functions
	CoupledConvectionReactionSub (const InputParameters &parameters)

const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialProperty (const std::string &name)

const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialPropertyByName (const std::string &name)

void	validateDerivativeMaterialPropertyBase (const std::string &base)

const MaterialPropertyName	derivativePropertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	derivativePropertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	derivativePropertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	derivativePropertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)

void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)

void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

Static Public Member Functions
static InputParameters	validParams ()

Protected Member Functions
virtual Real	computeQpResidual () override

virtual Real	computeQpJacobian () override

virtual Real	computeQpOffDiagJacobian (unsigned int jvar) override

Protected Attributes
const Real	_weight
	Weight of the equilibrium species concentration in the total primary species concentration.

const VariableValue &	_log_k
	Equilibrium constant for the equilibrium species in association form.

const Real	_sto_u
	Stoichiometric coefficient of the primary species.

const std::vector< Real >	_sto_v
	Stoichiometric coefficients of the coupled primary species.

const MaterialProperty< Real > &	_cond
	Hydraulic conductivity.

const RealVectorValue	_gravity
	Gravity.

const MaterialProperty< Real > &	_density
	Fluid density.

const VariableGradient &	_grad_p
	Pressure gradient.

const unsigned int	_pvar
	Pressure variable number.

const std::vector< unsigned int >	_vars
	Coupled primary species variable numbers.

const std::vector< const VariableValue * >	_vals
	Coupled primary species concentrations.

const std::vector< const VariableGradient * >	_grad_vals
	Coupled gradients of primary species concentrations.

const VariableValue &	_gamma_u
	Activity coefficient of primary species in the equilibrium species.

const std::vector< const VariableValue * >	_gamma_v
	Activity coefficients of coupled primary species in the equilibrium species.

const VariableValue &	_gamma_eq
	Activity coefficient of equilibrium species.

Private Member Functions
bool	haveMaterialProperty (const std::string &prop_name)

std::vector< VariableName >	buildVariableVector (const VariableName &c1, const VariableName &c2, const VariableName &c3)

void	validateCouplingHelper (const MaterialPropertyName &base, const std::vector< VariableName > &c, const System &system, std::vector< VariableName > &missing)

bool	isNotObjectVariable (const VariableName &name)

Private Attributes
FEProblemBase &	_dmi_fe_problem

Detailed Description

Convection of primary species in given equilibrium species.

Definition at line 18 of file CoupledConvectionReactionSub.h.

Constructor & Destructor Documentation

◆ CoupledConvectionReactionSub()

CoupledConvectionReactionSub::CoupledConvectionReactionSub ( const InputParameters & parameters )

Definition at line 40 of file CoupledConvectionReactionSub.C.

  : DerivativeMaterialInterface<Kernel>(parameters),
    _weight(getParam<Real>("weight")),
    _log_k(coupledValue("log_k")),
    _sto_u(getParam<Real>("sto_u")),
    _sto_v(getParam<std::vector<Real>>("sto_v")),
    _cond(getMaterialProperty<Real>("conductivity")),
    _gravity(getParam<RealVectorValue>("gravity")),
    _density(getDefaultMaterialProperty<Real>("density")),
    _grad_p(coupledGradient("p")),
    _pvar(coupled("p")),
    _vars(coupledIndices("v")),
    _vals(coupledValues("v")),
    _grad_vals(coupledGradients("v")),
    _gamma_u(coupledValue("gamma_u")),
    _gamma_v(isCoupled("gamma_v")
                 ? coupledValues("gamma_v") // have value
                 : std::vector<const VariableValue *>(coupledComponents("v"),
                                                      &coupledValue("gamma_v"))), // default
    _gamma_eq(coupledValue("gamma_eq"))
{
  const unsigned int n = coupledComponents("v");
 
  // Check that the correct number of coupled values have been provided
  if (_sto_v.size() != n)
    mooseError("The number of stoichiometric coefficients in sto_v is not equal to the number of "
               "coupled species in ",
               _name);
 
  if (isCoupled("gamma_v"))
    if (coupledComponents("gamma_v") != n)
      mooseError("The number of activity coefficients in gamma_v is not equal to the number of "
                 "coupled species in ",
                 _name);
}

Member Function Documentation

◆ computeQpJacobian()

Real CoupledConvectionReactionSub::computeQpJacobian ( )

overrideprotectedvirtual

Definition at line 106 of file CoupledConvectionReactionSub.C.

{
  RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
 
  RealGradient d_u_1 =
      _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_phi[_j][_qp];
  RealGradient d_u_2 = _phi[_j][_qp] * _sto_u * (_sto_u - 1.0) * _gamma_u[_qp] * _gamma_u[_qp] *
                       std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 2.0) * _grad_u[_qp];
 
  RealGradient d_var_sum(0.0, 0.0, 0.0);
  const Real d_v_u =
      _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _phi[_j][_qp];
 
  for (unsigned int i = 0; i < _vals.size(); ++i)
  {
    d_u_1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
    d_u_2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
 
    RealGradient d_var = d_v_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
                         std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
                         (*_grad_vals[i])[_qp];
    for (unsigned int j = 0; j < _vals.size(); ++j)
      if (j != i)
        d_var *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
 
    d_var_sum += d_var;
  }
 
  RealGradient d_u_j = d_u_1 + d_u_2;
  return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * darcy_vel * (d_u_j + d_var_sum) /
         _gamma_eq[_qp];
}

◆ computeQpOffDiagJacobian()

Real CoupledConvectionReactionSub::computeQpOffDiagJacobian ( unsigned int jvar )

overrideprotectedvirtual

Definition at line 140 of file CoupledConvectionReactionSub.C.

{
  if (jvar == _pvar)
  {
    RealVectorValue ddarcy_vel_dp = -_cond[_qp] * _grad_phi[_j][_qp];
 
    RealGradient d_u =
        _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
    RealGradient d_var_sum(0.0, 0.0, 0.0);
    const Real d_v_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
 
    for (unsigned int i = 0; i < _vals.size(); ++i)
    {
      d_u *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
 
      RealGradient d_var = d_v_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
                           std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
                           (*_grad_vals[i])[_qp];
      for (unsigned int j = 0; j < _vals.size(); ++j)
        if (j != i)
          d_var *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
 
      d_var_sum += d_var;
    }
    return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * ddarcy_vel_dp *
           (d_u + d_var_sum) / _gamma_eq[_qp];
  }
 
  if (_vals.size() == 0)
    return 0.0;
 
  RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
  RealGradient diff1 =
      _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
  for (unsigned int i = 0; i < _vals.size(); ++i)
  {
    if (jvar == _vars[i])
      diff1 *= _sto_v[i] * (*_gamma_v[i])[_qp] *
               std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
    else
      diff1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
  }
 
  Real val_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
  RealGradient diff2_1(1.0, 1.0, 1.0);
  RealGradient diff2_2(1.0, 1.0, 1.0);
  for (unsigned int i = 0; i < _vals.size(); ++i)
    if (jvar == _vars[i])
    {
      diff2_1 = _sto_v[i] * (*_gamma_v[i])[_qp] * (*_gamma_v[i])[_qp] * (_sto_v[i] - 1.0) *
                std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 2.0) * _phi[_j][_qp] *
                (*_grad_vals[i])[_qp];
      diff2_2 = _sto_v[i] * (*_gamma_v[i])[_qp] *
                std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
                _grad_phi[_j][_qp];
    }
 
  RealGradient diff2 = val_u * (diff2_1 + diff2_2);
  for (unsigned int i = 0; i < _vals.size(); ++i)
    if (jvar != _vars[i])
    {
      diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
    }
 
  RealGradient diff3;
  RealGradient diff3_sum(0.0, 0.0, 0.0);
  Real val_jvar = 0.0;
  unsigned int var = 0;
 
  for (unsigned int i = 0; i < _vals.size(); ++i)
    if (jvar == _vars[i])
    {
      var = i;
      val_jvar = val_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
                 std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) * _phi[_j][_qp];
    }
 
  for (unsigned int i = 0; i < _vals.size(); ++i)
    if (i != var)
    {
      diff3 = val_jvar * _sto_v[i] * (*_gamma_v[i])[_qp] *
              std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
              (*_grad_vals[i])[_qp];
 
      for (unsigned int j = 0; j < _vals.size(); ++j)
        if (j != var && j != i)
          diff3 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
 
      diff3_sum += diff3;
    }
 
  return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * darcy_vel *
         (diff1 + diff2 + diff3_sum) / _gamma_eq[_qp];
}

◆ computeQpResidual()

Real CoupledConvectionReactionSub::computeQpResidual ( )

overrideprotectedvirtual

Definition at line 77 of file CoupledConvectionReactionSub.C.

{
  RealVectorValue darcy_vel = -_cond[_qp] * (_grad_p[_qp] - _density[_qp] * _gravity);
  RealGradient d_u =
      _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_u[_qp];
  RealGradient d_var_sum(0.0, 0.0, 0.0);
  const Real d_v_u = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
 
  for (unsigned int i = 0; i < _vals.size(); ++i)
  {
    d_u *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
 
    RealGradient d_var = d_v_u * _sto_v[i] * (*_gamma_v[i])[_qp] *
                         std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i] - 1.0) *
                         (*_grad_vals[i])[_qp];
 
    for (unsigned int j = 0; j < _vals.size(); ++j)
      if (j != i)
        d_var *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
 
    d_var_sum += d_var;
  }
 
  mooseAssert(_gamma_eq[_qp] > 0.0, "Activity coefficient must be greater than zero");
  return _weight * std::pow(10.0, _log_k[_qp]) * _test[_i][_qp] * darcy_vel * (d_u + d_var_sum) /
         _gamma_eq[_qp];
}

◆ validParams()

InputParameters CoupledConvectionReactionSub::validParams ( )

static

Definition at line 15 of file CoupledConvectionReactionSub.C.

{
  InputParameters params = Kernel::validParams();
  params.addParam<Real>("weight", 1.0, "Weight of the equilibrium species");
  params.addCoupledVar("log_k", 0.0, "Equilibrium constant of dissociation equilibrium reaction");
  params.addParam<Real>("sto_u",
                        1.0,
                        "Stoichiometric coef of the primary species the kernel "
                        "operates on in the equilibrium reaction");
  params.addCoupledVar(
      "gamma_u", 1.0, "Activity coefficient of primary species that this kernel operates on");
  params.addParam<std::vector<Real>>(
      "sto_v",
      {},
      "The stoichiometric coefficients of coupled primary species in equilibrium reaction");
  params.addRequiredCoupledVar("p", "Pressure");
  params.addCoupledVar("v", "List of coupled primary species");
  params.addCoupledVar("gamma_v", 1.0, "Activity coefficients of coupled primary species");
  params.addCoupledVar("gamma_eq", 1.0, "Activity coefficient of this equilibrium species");
  RealVectorValue g(0, 0, 0);
  params.addParam<RealVectorValue>("gravity", g, "Gravity vector (default is (0, 0, 0))");
  params.addClassDescription("Convection of equilibrium species");
  return params;
}