Merge pull request #27737 from parmar0/materialInterfaceReaction

Add a material interface reaction

framework/doc/content/source/interfacekernels/ADMatInterfaceReaction.md

@@ -0,0 +1,41 @@
+# ADMatInterfaceReaction
+
+## 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}
+\textrm{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. `ADMatInterfaceReaction` object is used to impose this condition.
+
+[InterfaceDiffusion.md] is also used in this case to control flux at the interface.
+
+\begin{equation}
+D_1 \frac {\partial C_1} {\partial n} = D_2 \frac {\partial C_2} {\partial n}
+\end{equation}
+
+However, the flux is not [well-defined](https://en.wikipedia.org/wiki/Well-defined_expression) across the interface. The `ADMatInterfaceReaction` interfacekernel applies a condition to constrain the potential discontinuity across the interface.
+
+Both kernels at the interface work together to give full mathematical and physical meaning of the problem. Together, the implicit equations represented by `ADMatInterfaceReaction` and [InterfaceDiffusion.md] combine to provide the following relationship at the interface.
+
+\begin{equation}
+D_1 \frac {\partial C_1} {\partial n} = D_2 \frac {\partial C_2} {\partial n} + k_f C_1 - k_b C_2
+\end{equation}
+
+## Example Input Syntax
+
+!listing test/tests/interfacekernels/1d_interface/ADMatreaction_1D_steady.i start=[interface_reaction] end=[] include-end=true
+
+!listing test/tests/interfacekernels/1d_interface/ADMatreaction_1D_transient.i start=[interface_reaction] end=[] include-end=true
+
+!syntax parameters /InterfaceKernels/InterfaceReaction
+
+!syntax children /InterfaceKernels/ADMatInterfaceReaction

framework/include/interfacekernels/ADMatInterfaceReaction.h

@@ -0,0 +1,33 @@
+//* 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 "ADInterfaceKernel.h"
+
+/**
+ * Implements a reaction to establish ReactionRate=k_f*u-k_b*v
+ * at interface.
+ */
+class ADMatInterfaceReaction : public ADInterfaceKernel
+{
+public:
+  static InputParameters validParams();
+
+  ADMatInterfaceReaction(const InputParameters & parameters);
+
+protected:
+  virtual ADReal computeQpResidual(Moose::DGResidualType type) override;
+
+  /// Forward reaction rate coefficient
+  const ADMaterialProperty<Real> & _kf;
+
+  /// Backward reaction rate coefficient
+  const ADMaterialProperty<Real> & _kb;
+};

framework/src/interfacekernels/ADMatInterfaceReaction.C

@@ -0,0 +1,54 @@
+//* 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 "ADMatInterfaceReaction.h"
+
+registerMooseObject("MooseApp", ADMatInterfaceReaction);
+
+InputParameters
+ADMatInterfaceReaction::validParams()
+
+{
+  InputParameters params = ADInterfaceKernel::validParams();
+  params.addParam<MaterialPropertyName>("forward_rate", "kf", "Forward reaction rate coefficient.");
+  params.addParam<MaterialPropertyName>(
+      "backward_rate", "kb", "Backward reaction rate coefficient.");
+  params.addClassDescription("Implements a reaction to establish ReactionRate=k_f*u-k_b*v "
+                             "at interface.");
+  return params;
+}
+
+ADMatInterfaceReaction::ADMatInterfaceReaction(const InputParameters & parameters)
+  : ADInterfaceKernel(parameters),
+    _kf(getADMaterialProperty<Real>("forward_rate")),
+    _kb(getNeighborADMaterialProperty<Real>("backward_rate"))
+{
+}
+
+ADReal
+ADMatInterfaceReaction::computeQpResidual(Moose::DGResidualType type)
+{
+  ADReal r = 0;
+  switch (type)
+  {
+    // Move all the terms to the LHS to get residual, for primary domain
+    // Residual = kf*u - kb*v = kf*u - kb*v
+    // Weak form for primary domain is: (test, kf*u - kb*v)
+    case Moose::Element:
+      r = _test[_i][_qp] * (_kf[_qp] * _u[_qp] - _kb[_qp] * _neighbor_value[_qp]);
+      break;
+
+    // Similarly, weak form for secondary domain is: -(test, kf*u - kb*v),
+    // flip the sign because the direction is opposite.
+    case Moose::Neighbor:
+      r = -_test_neighbor[_i][_qp] * (_kf[_qp] * _u[_qp] - _kb[_qp] * _neighbor_value[_qp]);
+      break;
+  }
+  return r;
+}

test/tests/interfacekernels/1d_interface/ADMatreaction_1D_steady.i

@@ -0,0 +1,170 @@
+# Steady-state test for the ADMatInterfaceReaction 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 = 0
+#
+# This results in the interfacial conditions
+#
+# -D1 du = -D2 dv
+#   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]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 1
+    nx = 10
+    xmax = 2
+  []
+  [subdomain1]
+    input = gen
+    type = SubdomainBoundingBoxGenerator
+    bottom_left = '1.0 0 0'
+    block_id = 1
+    top_right = '2.0 1.0 0'
+  []
+  [interface]
+    type = SideSetsBetweenSubdomainsGenerator
+    input = 'subdomain1'
+    primary_block = '0'
+    paired_block = '1'
+    new_boundary = 'primary0_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'
+    diffusivity = D
+  []
+  [diff_v]
+    type = MatDiffusion
+    variable = v
+    block = '1'
+    diffusivity = D
+  []
+[]
+
+[InterfaceKernels]
+  [interface]
+    type = InterfaceDiffusion
+    variable = u
+    neighbor_var = 'v'
+    boundary = 'primary0_interface'
+    D = D
+    D_neighbor = D
+  []
+  [interface_reaction]
+    type = ADMatInterfaceReaction
+    variable = u
+    neighbor_var = 'v'
+    boundary = 'primary0_interface'
+    forward_rate = forward_rate
+    backward_rate = backward_rate
+  []
+[]
+
+[BCs]
+  [left]
+    type = ADDirichletBC
+    variable = u
+    boundary = 'left'
+    value = 1
+  []
+  [right]
+    type = ADDirichletBC
+    variable = v
+    boundary = 'right'
+    value = 0
+  []
+[]
+
+[Materials]
+  [block0]
+    type = 'ADGenericConstantMaterial'
+    block = '0'
+    prop_names = 'forward_rate backward_rate'
+    prop_values = '1.0 2.0'
+  []
+  [block01]
+    type = 'GenericConstantMaterial'
+    block = '0'
+    prop_names = 'D'
+    prop_values = '4'
+  []
+  [block1]
+    type = 'ADGenericConstantMaterial'
+    block = '1'
+    prop_names = 'forward_rate backward_rate'
+    prop_values = '1.0 2.0'
+  []
+  [block11]
+    type = 'GenericConstantMaterial'
+    block = '1'
+    prop_names = 'D'
+    prop_values = '2'
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+  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'
+  []
+[]