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

#include <CoupledConvectionReactionSub.h>

Inheritance diagram for CoupledConvectionReactionSub:

Public Member Functions
	CoupledConvectionReactionSub (const InputParameters &parameters)

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

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

const Real	_sto_u
	Stoichiometric coefficient of the primary species. More...

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

const MaterialProperty< Real > &	_cond
	Hydraulic conductivity. More...

const RealVectorValue	_gravity
	Gravity. More...

const MaterialProperty< Real > &	_density
	Fluid density. More...

const VariableGradient &	_grad_p
	Pressure gradient. More...

const unsigned int	_pvar
	Pressure variable number. More...

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

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

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

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

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

const VariableValue &	_gamma_eq
	Activity coefficient of equilibrium species. More...

Detailed Description

Convection of primary species in given equilibrium species.

Definition at line 23 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")),
     _gamma_u(coupledValue("gamma_u")),
     _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);
   _vars.resize(n);
   _vals.resize(n);
   _grad_vals.resize(n);
   _gamma_v.resize(n);
  
   for (unsigned int i = 0; i < n; ++i)
   {
     _vars[i] = coupled("v", i);
     _vals[i] = &coupledValue("v", i);
     _grad_vals[i] = &coupledGradient("v", i);
     // If gamma_v has been supplied, use those values, but if not, use the default value
     _gamma_v[i] = (isCoupled("gamma_v") ? &coupledValue("gamma_v", i) : &coupledValue("gamma_v"));
   }
 }

Member Function Documentation

◆ computeQpJacobian()

Real CoupledConvectionReactionSub::computeQpJacobian ( )

overrideprotectedvirtual

Definition at line 112 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 146 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 83 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];
 }