diff --git a/include/interfacekernels/InterfaceSorption.h b/include/interfacekernels/InterfaceSorption.h index 898b34102..c6480653a 100644 --- a/include/interfacekernels/InterfaceSorption.h +++ b/include/interfacekernels/InterfaceSorption.h @@ -62,6 +62,11 @@ class InterfaceSorptionTempl : public InterfaceKernelParent const GenericMaterialProperty * _diffusivity_neighbor; ///@} + + /// Defines whether a factor of RT is needeed + const bool _neighbour_is_gas; + + ///@{Retrieve parent members using InterfaceKernelParent::_u; using InterfaceKernelParent::_qp; diff --git a/src/interfacekernels/InterfaceSorption.C b/src/interfacekernels/InterfaceSorption.C index 735ce56bf..070fbb62c 100644 --- a/src/interfacekernels/InterfaceSorption.C +++ b/src/interfacekernels/InterfaceSorption.C @@ -51,8 +51,13 @@ InterfaceSorptionTempl::validParams() false, "Whether to use the penalty formulation to enforce the flux balance"); params.addParam("flux_penalty", 1.0, "Penalty associated with the flux balance"); - params.addParam("diffusivity", "diffusivity", "The diffusion coefficient"); + + + params.addParam("neighbour_is_gas", + true, + "When the neighbour concentration is a gas a factor of RT is applied."); + return params; } @@ -61,23 +66,19 @@ InterfaceSorptionTempl::InterfaceSorptionTempl(const InputParameters & pa : InterfaceKernelParent(parameters), _K0(this->template getParam("K0")), _Ea(this->template getParam("Ea")), - _n_sorption(this->template getParam("n_sorption")), - _unit_scale(this->template getParam("unit_scale")), _unit_scale_neighbor(this->template getParam("unit_scale_neighbor")), - _T(this->template coupledGenericValue("temperature")), - _sorption_penalty(this->template getParam("sorption_penalty")), _use_flux_penalty(this->template getParam("use_flux_penalty")), _flux_penalty(_use_flux_penalty ? this->template getParam("flux_penalty") : 0.0), - _diffusivity(this->template getGenericMaterialProperty("diffusivity")), _diffusivity_neighbor( - _use_flux_penalty - ? &this->template getGenericNeighborMaterialProperty("diffusivity") - : nullptr) + _use_flux_penalty + ? &this->template getGenericNeighborMaterialProperty("diffusivity") + : nullptr), + _neighbour_is_gas(this->template getParam("neighbour_is_gas")) { if (!_use_flux_penalty && parameters.isParamSetByUser("flux_penalty")) this->template paramError( @@ -101,13 +102,16 @@ InterfaceSorptionTempl::computeQpResidual(Moose::DGResidualType type) // The unit scaling affects the residual of the sorption equation, but not the flux equation. // Otherwise mass is not conserved. + // For a gas: p = RT c + auto RT_scalefactor = _neighbour_is_gas ? (R * temperature_limited) : 1.0; + switch (type) { case Moose::Element: { residual = _test[_i][_qp] * _sorption_penalty * (u - _K0 * std::exp(-_Ea / R / temperature_limited) * - std::pow(u_neighbor * R * temperature_limited, _n_sorption)); + std::pow(u_neighbor * RT_scalefactor, _n_sorption)); break; } @@ -140,6 +144,8 @@ InterfaceSorptionTempl::computeQpJacobian(Moose::DGJacobianType type) // The unit scaling affects the jacobian of the sorption equation, but not the flux equation. // Otherwise mass is not conserved. + auto RT_scalefactor = _neighbour_is_gas ? (R * temperature_limited) : 1.0; + switch (type) { case Moose::ElementElement: @@ -150,7 +156,7 @@ InterfaceSorptionTempl::computeQpJacobian(Moose::DGJacobianType type) jacobian = -_test[_i][_qp] * _sorption_penalty * _unit_scale_neighbor / _unit_scale * _phi_neighbor[_j][_qp] * (_K0 * std::exp(-_Ea / R / temperature_limited) * - std::pow(R * temperature_limited, _n_sorption) * _n_sorption * + std::pow(RT_scalefactor, _n_sorption) * _n_sorption * std::pow(u_neighbor, _n_sorption - 1.)); break;