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

#include <CoupledDiffusionReactionSub.h>

Inheritance diagram for CoupledDiffusionReactionSub:

Public Member Functions
	CoupledDiffusionReactionSub (const InputParameters &parameters)

Protected Member Functions
virtual Real	computeQpResidual () override

virtual Real	computeQpJacobian () override

virtual Real	computeQpOffDiagJacobian (unsigned int jvar) override

Private Attributes
const MaterialProperty< Real > &	_diffusivity
	Material property of dispersion-diffusion coefficient. More...

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

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

Diffusion of primary species in given equilibrium species.

Definition at line 22 of file CoupledDiffusionReactionSub.h.

Constructor & Destructor Documentation

◆ CoupledDiffusionReactionSub()

CoupledDiffusionReactionSub::CoupledDiffusionReactionSub ( const InputParameters & parameters )

Definition at line 40 of file CoupledDiffusionReactionSub.C.

   : Kernel(parameters),
     _diffusivity(getMaterialProperty<Real>("diffusivity")),
     _weight(getParam<Real>("weight")),
     _log_k(coupledValue("log_k")),
     _sto_u(getParam<Real>("sto_u")),
     _sto_v(getParam<std::vector<Real>>("sto_v")),
     _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 CoupledDiffusionReactionSub::computeQpJacobian ( )

overrideprotectedvirtual

Definition at line 109 of file CoupledDiffusionReactionSub.C.

 {
   RealGradient diff1_1 =
       _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _grad_phi[_j][_qp];
   RealGradient diff1_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];
   for (unsigned int i = 0; i < _vals.size(); ++i)
   {
     diff1_1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
     diff1_2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
   }
  
   RealGradient diff1 = diff1_1 + diff1_2;
   Real d_val =
       _sto_u * _gamma_u[_qp] * std::pow(_gamma_u[_qp] * _u[_qp], _sto_u - 1.0) * _phi[_j][_qp];
   RealGradient diff2_sum(0.0, 0.0, 0.0);
   for (unsigned int i = 0; i < _vals.size(); ++i)
   {
     RealGradient diff2 = d_val * _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)
         diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
  
     diff2_sum += diff2;
   }
  
   return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
          (diff1 + diff2_sum) / _gamma_eq[_qp];
 }

◆ computeQpOffDiagJacobian()

Real CoupledDiffusionReactionSub::computeQpOffDiagJacobian ( unsigned int jvar )

overrideprotectedvirtual

Definition at line 142 of file CoupledDiffusionReactionSub.C.

 {
   // If no coupled species, return 0
   if (_vals.size() == 0)
     return 0.0;
  
   // If jvar is not one of the coupled species, return 0
   if (std::find(_vars.begin(), _vars.end(), jvar) == _vars.end())
     return 0.0;
  
   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] * (_sto_v[i] - 1.0) * (*_gamma_v[i])[_qp] * (*_gamma_v[i])[_qp] *
                 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]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
          (diff1 + diff2 + diff3_sum) / _gamma_eq[_qp];
 }

◆ computeQpResidual()

Real CoupledDiffusionReactionSub::computeQpResidual ( )

overrideprotectedvirtual

Definition at line 80 of file CoupledDiffusionReactionSub.C.

 {
   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)
     diff1 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[i])[_qp], _sto_v[i]);
  
   RealGradient diff2_sum(0.0, 0.0, 0.0);
   const Real d_val = std::pow(_gamma_u[_qp] * _u[_qp], _sto_u);
   for (unsigned int i = 0; i < _vals.size(); ++i)
   {
     RealGradient diff2 = d_val * _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)
         diff2 *= std::pow((*_gamma_v[i])[_qp] * (*_vals[j])[_qp], _sto_v[j]);
  
     diff2_sum += diff2;
   }
  
   mooseAssert(_gamma_eq[_qp] > 0.0, "Activity coefficient must be greater than zero");
   return _weight * std::pow(10.0, _log_k[_qp]) * _diffusivity[_qp] * _grad_test[_i][_qp] *
          (diff1 + diff2_sum) / _gamma_eq[_qp];
   ;
 }