Skip to content
Permalink
Browse files

Merge pull request #13366 from wesleyzzz/Test

Add Interface Reaction Kernel and examples
  • Loading branch information...
lindsayad committed Jun 3, 2019
2 parents 62ac1db + ea31281 commit 3ed1345228318613b52512948e54caa7c3ce09c7
@@ -0,0 +1,42 @@
# InterfaceReaction

## Description

Specie M transports between two domains (domain 1 and domain 2), at the interface consider the following reaction is taking place:

\begin{equation}
M(1)\xrightleftharpoons[k_b]{k_f}M(2)
\end{equation}

With the first order reaction rate assuming a quasi-steady-state

\begin{equation}
Reaction Rate = \frac {\partial C_1} {\partial t} = k_f C_1 - k_b C_2 \approx 0
\end{equation}

where $C_1$ is the specie concentration in domain 1, $C_2$ is the specie concentration in domain 2, $k_f$ is the forward reaction coefficient, and $k_b$ is the backward reaction coefficient. `InterfaceReaction` object is used to impose this condition. Associated kernel is:

[/InterfaceReaction.C]

[/InterfaceReaction.h]

In addition, fluxes are matched from both domains, this could be achieved by [`InterfaceDiffusion`](/InterfaceKernels/index.md).

Both kernels at the interface work together to give full mathematical and physical meaning of the problem.

Two examples (steady-state and transient-state) are shown in the MOOSE test directory,

[1d_interface/reaction_1D_steady.i]

[1d_interface/reaction_1D_transient.i]


## Example Input Syntax

!listing test/tests/interfacekernels/1d_interface/reaction_1D_steady.i start=[./interface_reaction] end=[../] include-end=true

!syntax parameters /InterfaceKernels/InterfaceReaction

!syntax inputs /InterfaceKernels/InterfaceReaction

!syntax children /InterfaceKernels/InterfaceReaction
@@ -0,0 +1,38 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "InterfaceKernel.h"

// Forward Declarations
class InterfaceReaction;

template <>
InputParameters validParams<InterfaceReaction>();

/**
* Implements a reaction to establish ReactionRate=k_f*u-k_b*v
* at interface.
*/
class InterfaceReaction : public InterfaceKernel
{
public:
InterfaceReaction(const InputParameters & parameters);

protected:
virtual Real computeQpResidual(Moose::DGResidualType type) override;
virtual Real computeQpJacobian(Moose::DGJacobianType type) override;

/// Forward reaction rate coefficient
Real _kf;

/// Backward reaction rate coefficient
Real _kb;
};
@@ -0,0 +1,73 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "InterfaceReaction.h"

registerMooseObject("MooseApp", InterfaceReaction);

template <>
InputParameters
validParams<InterfaceReaction>()
{
InputParameters params = validParams<InterfaceKernel>();
params.addRequiredParam<Real>("kf", "Forward reaction rate coefficient.");
params.addRequiredParam<Real>("kb", "Backward reaction rate coefficient.");
params.addClassDescription("Implements a reaction to establish ReactionRate=k_f*u-k_b*v "
"at interface.");
return params;
}

InterfaceReaction::InterfaceReaction(const InputParameters & parameters)
: InterfaceKernel(parameters), _kf(getParam<Real>("kf")), _kb(getParam<Real>("kb"))
{
}

Real
InterfaceReaction::computeQpResidual(Moose::DGResidualType type)
{
Real r = 0;
switch (type)
{
// Move all the terms to the LHS to get residual, for master domain
// Residual = kf*u - kb*v = kf*u - kb*v
// Weak form for master domain is: (test, kf*u - kb*v)
case Moose::Element:
r = _test[_i][_qp] * (_kf * _u[_qp] - _kb * _neighbor_value[_qp]);
break;

// Similarly, weak form for slave domain is: -(test, kf*u - kb*v),
// flip the sign because the direction is opposite.
case Moose::Neighbor:
r = -_test_neighbor[_i][_qp] * (_kf * _u[_qp] - _kb * _neighbor_value[_qp]);
break;
}
return r;
}

Real
InterfaceReaction::computeQpJacobian(Moose::DGJacobianType type)
{
Real jac = 0;
switch (type)
{
case Moose::ElementElement:
jac = _test[_i][_qp] * _kf * _phi[_j][_qp];
break;
case Moose::NeighborNeighbor:
jac = -_test_neighbor[_i][_qp] * -_kb * _phi_neighbor[_j][_qp];
break;
case Moose::NeighborElement:
jac = -_test_neighbor[_i][_qp] * _kf * _phi[_j][_qp];
break;
case Moose::ElementNeighbor:
jac = _test[_i][_qp] * -_kb * _phi_neighbor[_j][_qp];
break;
}
return jac;
}
@@ -0,0 +1,2 @@
time,elemental_error_u,elemental_error_v
2,7.0216669371534e-17,3.4152471690487e-17
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,160 @@
# Steady-state test for the InterfaceReaction kernel.
#
# Specie M transport from domain 1 (0<=x<=1) to domain 2 (1<x<=2),
# u and v are concentrations in domain 1 and domain 2.
#
# Diffusion in both domains can be described by Ficks law and diffusion
# kernel is applied.
#
# Specie M has different diffusity in different domains, here set as D1=4, D2=2.
#
# Dirichlet boundary conditions are applied, i.e., u(0)=1, v(2)=0
#
# At the interface consider the following
#
# (a) Fluxes are matched from both domains (InterfaceDiffusion kernel)
#
# (b) First-order reaction is R = kf*u - kb*v
#
# Analytical solution is
# u = -0.2*u+1, 0<=u<=1
# v = -0.4*v+0.8, 1<v<=2

[Mesh]
type = GeneratedMesh
dim = 1
nx = 10
xmax = 2
[]

[MeshModifiers]
[./subdomain1]
type = SubdomainBoundingBox
bottom_left = '1.0 0 0'
block_id = 1
top_right = '2.0 1.0 0'
[../]
[./interface]
type = SideSetsBetweenSubdomains
depends_on = 'subdomain1'
master_block = '0'
paired_block = '1'
new_boundary = 'master0_interface'
[../]
[]

[Variables]
[./u]
order = FIRST
family = LAGRANGE
block = '0'
[../]
[./v]
order = FIRST
family = LAGRANGE
block = '1'
[../]
[]

[Kernels]
[./diff_u]
type = MatDiffusion
variable = u
block = '0'
D_name = D
[../]
[./diff_v]
type = MatDiffusion
variable = v
block = '1'
D_name = D
[../]
[]

[InterfaceKernels]
[./interface]
type = InterfaceDiffusion
variable = u
neighbor_var = 'v'
boundary = 'master0_interface'
D = D
D_neighbor = D
[../]
[./interface_reaction]
type = InterfaceReaction
variable = u
neighbor_var = 'v'
boundary = 'master0_interface'
kf = 1 # Forward reaction rate coefficient
kb = 2 # Backward reaction rate coefficient
[../]
[]

[BCs]
[./left]
type = DirichletBC
variable = u
boundary = 'left'
value = 1
[../]
[./right]
type = DirichletBC
variable = v
boundary = 'right'
value = 0
[../]
[]

[Materials]
[./block0]
type = GenericConstantMaterial
block = '0'
prop_names = 'D'
prop_values = '4'
[../]
[./block1]
type = GenericConstantMaterial
block = '1'
prop_names = 'D'
prop_values = '2'
[../]
[]

[Preconditioning]
[./smp]
type = SMP
full = true
[../]
[]

[Executioner]
type = Steady
solve_type = PJFNK
nl_rel_tol = 1e-10
[]

[Outputs]
print_linear_residuals = true
execute_on = 'FINAL'
exodus = true
csv = true
[]

[Debug]
show_var_residual_norms = true
[]

[Postprocessors]
[./elemental_error_u]
type = ElementL2Error
function = -0.2*x+1
variable = 'u'
block = '0'
[../]
[./elemental_error_v]
type = ElementL2Error
function = -0.4*x+0.8
variable = 'v'
block = '1'
[../]
[]

0 comments on commit 3ed1345

Please sign in to comment.
You can’t perform that action at this time.