Repository: PMEAL/OpenPNM
Pull request #1492 (2020-05-28 15:44:07)

Pore-scale models for Butler-Volmer kinetics
This PR adds two pore-scale models for Butler-Volmer kinetics; one for concentration to be used in mass transfer algorithms, and the other for voltage to be used in electron conduction algorithms. The two pore-scale models are complementary, meaning that they must be used together.

Here's an example script on how to use them:
```python
# Define parameters for Butler-Volmer pore-scale model
BV_params = {
    "z": 4,
    "j0": 1e-3,
    "c_ref": 1000,
    "alpha_anode": 0.4,
    "alpha_cathode": 0.6,
    "reaction_order": 1
}
BV_c = op.models.physics.generic_source_term.butler_volmer_conc
BV_v = op.models.physics.generic_source_term.butler_volmer_voltage

# Add source term for O2 consumption
phys.add_model(propname="pore.rxn_fd", model=BV_c,
               X="pore.concentration",
               reaction_area="pore.reaction_area",
               solid_voltage="pore.solid_voltage",
               electrolyte_voltage="pore.voltage",
               open_circuit_voltage="pore.open_circuit_voltage",
               regen_mode="deferred", **BV_params)

# Add source term for H+ consumption
phys.add_model(propname="pore.rxn_oc", model=BV_v,
               X="pore.voltage",
               reaction_area="pore.reaction_area",
               solid_voltage="pore.solid_voltage",
               electrolyte_concentration="pore.concentration",
               open_circuit_voltage="pore.open_circuit_voltage",
               regen_mode="deferred", **BV_params)
```
Alternatively, you could explicitly pass parameters inside `BC_params` dictionary as follows, but I find the above approach more appealing and convenient.
```python
phys.add_model(propname="pore.rxn_oc", model=BV_v,
               X="pore.voltage",
               z=4, j0=1e-3, c_ref=1000, alpha_anode=0.4, alpha_cathode=0.6,
               reaction_order=1,
               reaction_area="pore.reaction_area",
               solid_voltage="pore.solid_voltage",
               electrolyte_concentration="pore.concentration",
               open_circuit_voltage="pore.open_circuit_voltage",
               regen_mode="deferred")
```


Patch:
diff --git a/example.py b/example.py
index a3175707c..ef48665ca 100644
--- a/example.py
+++ b/example.py
@@ -27,7 +27,7 @@ water.update(perm.results())
 # print(perm.calc_effective_permeability())
 
 # Add reaction term to phys_air
-mod = op.models.physics.generic_source_term.standard_kinetics
+mod = op.models.physics.source_terms.standard_kinetics
 phys_air['pore.n'] = 2
 phys_air['pore.A'] = -1e-5
 phys_air.add_model(propname='pore.2nd_order_rxn', model=mod,
diff --git a/openpnm/models/physics/__init__.py b/openpnm/models/physics/__init__.py
index ee0f0eea7..f02f9b36e 100644
--- a/openpnm/models/physics/__init__.py
+++ b/openpnm/models/physics/__init__.py
@@ -17,7 +17,8 @@ from . import thermal_conductance
 from . import hydraulic_conductance
 from . import dispersive_conductance
 from . import multiphase
-from . import generic_source_term
+from . import source_terms
+from . import source_terms as generic_source_term
 from . import flow_shape_factors
 from . import poisson_shape_factors
 from . import meniscus
diff --git a/openpnm/models/physics/generic_source_term.py b/openpnm/models/physics/source_terms.py
similarity index 67%
rename from openpnm/models/physics/generic_source_term.py
rename to openpnm/models/physics/source_terms.py
index b04655580..87c7a6a85 100644
--- a/openpnm/models/physics/generic_source_term.py
+++ b/openpnm/models/physics/source_terms.py
@@ -2,6 +2,7 @@ r"""
 Pore-scale models to be used as source term (e.g. reaction) in algorithms.
 """
 import numpy as _np
+import scipy as _sp
 
 __all__ = [
     "charge_conservation",
@@ -18,7 +19,9 @@ __all__ = [
     "natural_exponential_sym",
     "logarithm_sym",
     "natural_logarithm_sym",
-    "general_symbolic"
+    "general_symbolic",
+    "butler_volmer_conc",
+    "butler_volmer_voltage"
 ]
 
 
@@ -102,18 +105,25 @@ def charge_conservation(target, phase, p_alg, e_alg, assumption):
 def standard_kinetics(target, X, prefactor, exponent):
     r"""
     Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
+    function at the given value of ``X``:
 
         .. math::
             r = A X^b
 
     Parameters
     ----------
-    ``prefactor, exponent`` : {string}
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
+
+    X : string
+        The dictionary key on the target object containing the the quantity
+        of interest
+
+    prefactor, exponent : {string}
         The dictionary keys on the target object containing the coefficients
         values to be used in the source term model
 
-    ``quantity`` : {string}
+    quantity : {string}
         The dictionary key on the target object containing the the quantity
         of interest
 
@@ -146,6 +156,10 @@ def standard_kinetics(target, X, prefactor, exponent):
 
 
 def _parse_args(target, key, default):
+    r"""
+    Parses the passed ``key``: returns ``default`` if ``key`` is empty,
+    otherwise fetches and returns the values from ``target``.
+    """
     if key == '':
         val = default
     else:
@@ -156,21 +170,24 @@ def _parse_args(target, key, default):
 def linear(target, X, A1='', A2=''):
     r"""
     Calculates the rate, as well as slope and intercept of the following
-    function at the given value of `X`:
+    function at the given value of ``X``:
 
         .. math::
-            r = A_{1}   X  +  A_{2}
+            r = A_{1} X + A_{2}
 
     Parameters
     ----------
-    A1 -> A2 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A2 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -185,7 +202,7 @@ def linear(target, X, A1='', A2=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     A = _parse_args(target=target, key=A1, default=0.0)
@@ -209,14 +226,17 @@ def power_law(target, X, A1='', A2='', A3=''):
 
     Parameters
     ----------
-    A1 -> A3 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A3 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -231,7 +251,7 @@ def power_law(target, X, A1='', A2='', A3=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     A = _parse_args(target=target, key=A1, default=0.0)
@@ -256,14 +276,17 @@ def exponential(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
 
     Parameters
     ----------
-    A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A6 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -278,7 +301,7 @@ def exponential(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     A = _parse_args(target=target, key=A1, default=0.0)
@@ -302,18 +325,21 @@ def natural_exponential(target, X, A1='', A2='', A3='', A4='', A5=''):
     function at the given value of `X`:
 
         .. math::
-            r =   A_{1} exp( A_{2}  x^{ A_{3} } + A_{4} )+ A_{5}
+            r = A_{1} exp( A_{2}  x^{ A_{3} } + A_{4} )+ A_{5}
 
     Parameters
     ----------
-    A1 -> A5 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A5 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -328,7 +354,7 @@ def natural_exponential(target, X, A1='', A2='', A3='', A4='', A5=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     A = _parse_args(target=target, key=A1, default=0.0)
@@ -351,18 +377,21 @@ def logarithm(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     function at the given value of `X`:
 
         .. math::
-            r =  A_{1}   Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
+            r =  A_{1} Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
 
     Parameters
     ----------
-    A1 -> A6 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A6 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -377,7 +406,7 @@ def logarithm(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     A = _parse_args(target=target, key=A1, default=0.0)
@@ -402,18 +431,21 @@ def natural_logarithm(target, X, A1='', A2='', A3='', A4='', A5=''):
     function at the given value of `X`:
 
         .. math::
-            r =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
+            r =   A_{1} Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
 
     Parameters
     ----------
-    A1 -> A5 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A5 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -428,7 +460,7 @@ def natural_logarithm(target, X, A1='', A2='', A3='', A4='', A5=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     A = _parse_args(target=target, key=A1, default=0.0)
@@ -466,18 +498,21 @@ def linear_sym(target, X, A1='', A2=''):
     function at the given value of *x*:
 
         .. math::
-            r = A_{1}   x  +  A_{2}
+            r = A_{1} x + A_{2}
 
     Parameters
     ----------
-    A1 -> A2 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
         The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A2 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -492,7 +527,7 @@ def linear_sym(target, X, A1='', A2=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     from sympy import symbols
@@ -519,18 +554,21 @@ def power_law_sym(target, X, A1='', A2='', A3=''):
     function at the given value of *x*:
 
         .. math::
-            r = A_{1}   x^{A_{2}}  +  A_{3}
+            r = A_{1} x^{A_{2}} + A_{3}
 
     Parameters
     ----------
-    A1 -> A3 : string
-        The dictionary keys on the target object containing the coefficients
-        values to be used in the source term model
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
 
     X : string
-        The dictionary key on the target objecxt containing the the quantity
+        The dictionary key on the target object containing the the quantity
         of interest
 
+    A1 -> A3 : string
+        The dictionary keys on the target object containing the coefficients
+        values to be used in the source term model
+
     Returns
     -------
     A dictionary containing the following three items:
@@ -545,7 +583,7 @@ def power_law_sym(target, X, A1='', A2='', A3=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     from sympy import symbols
@@ -573,18 +611,21 @@ def exponential_sym(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     function at the given value of *x*:
 
         .. math::
-            r =  A_{1} A_{2}^{( A_{3} x^{ A_{4} } + A_{5})} + A_{6}
+            r = A_{1} A_{2}^{( A_{3} x^{ A_{4} } + A_{5})} + A_{6}
 
     Parameters
     ----------
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
+
+    X : string
+        The dictionary key on the target object containing the the quantity
+        of interest
+
     A1 -> A6 : string
         The dictionary keys on the target object containing the coefficients
         values to be used in the source term model
 
-    X : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
     Returns
     -------
     A dictionary containing the following three items:
@@ -599,7 +640,7 @@ def exponential_sym(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     from sympy import symbols
@@ -630,18 +671,21 @@ def natural_exponential_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
     function at the given value of *x*:
 
         .. math::
-            r =   A_{1} exp( A_{2}  x^{ A_{3} } + A_{4} )+ A_{5}
+            r = A_{1} exp( A_{2} x^{ A_{3} } + A_{4} )+ A_{5}
 
     Parameters
     ----------
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
+
+    X : string
+        The dictionary key on the target object containing the the quantity
+        of interest
+
     A1 -> A6 : string
         The dictionary keys on the target object containing the coefficients
         values to be used in the source term model
 
-    X : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
     Returns
     -------
     A dictionary containing the following three items:
@@ -656,7 +700,7 @@ def natural_exponential_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     from sympy import symbols, exp
@@ -686,18 +730,21 @@ def logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     function at the given value of *x*:
 
         .. math::
-            r =  A_{1}   Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
+            r = A_{1} Log_{ A_{2} }( A_{3} x^{ A_{4} }+ A_{5})+ A_{6}
 
     Parameters
     ----------
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
+
+    X : string
+        The dictionary key on the target object containing the the quantity
+        of interest
+
     A1 -> A6 : string
         The dictionary keys on the target object containing the coefficients
         values to be used in the source term model
 
-    x : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
     Returns
     -------
     A dictionary containing the following three items:
@@ -712,7 +759,7 @@ def logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5='', A6=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     from sympy import symbols, log
@@ -743,18 +790,21 @@ def natural_logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
     function at the given value of *x*:
 
         .. math::
-            rate =   A_{1}  Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
+            rate = A_{1} Ln( A_{2} x^{ A_{3} }+ A_{4})+ A_{5}
 
     Parameters
     ----------
+    target : OpenPNM object
+        The OpenPNM object where the result will be applied.
+
+    X : string
+        The dictionary key on the target object containing the the quantity
+        of interest
+
     A1 -> A5 : string
         The dictionary keys on the target object containing the coefficients
         values to be used in the source term model
 
-    X : string or float/int or array/list
-        The dictionary key on the target objecxt containing the the quantity
-        of interest
-
     Returns
     -------
     A dictionary containing the following three items:
@@ -769,7 +819,7 @@ def natural_logarithm_sym(target, X, A1='', A2='', A3='', A4='', A5=''):
     current value of X as follow:
 
         .. math::
-            rate = S_{1}   X  +  S_{2}
+            rate = S_{1} X + S_{2}
 
     """
     from sympy import symbols, ln
@@ -856,3 +906,243 @@ def general_symbolic(target, eqn=None, arg_map=None):
     s2_val = s2(*data.values())
     values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
     return values
+
+
+def butler_volmer_conc(
+    target, X, z, j0, c_ref, alpha_anode, alpha_cathode,
+    reaction_order=1,
+    temperature="pore.temperature",
+    reaction_area="pore.reaction_area",
+    solid_voltage="pore.solid_voltage",
+    electrolyte_voltage="pore.electrolyte_voltage",
+    open_circuit_voltage="pore.open_circuit_voltage",
+):
+    r"""
+    Calculates the rate, slope and intercept of the Butler-Volmer kinetic
+    model based on **concentration** to be used in mass transfer
+    algorithms.
+
+        .. math::
+            r = j_0 A_{rxn} (\frac{ X }{ c_{ref} }) ^ {\nu}
+            \Big(
+                \exp(  \frac{\alpha_a z F}{RT} \eta )
+              - \exp( -\frac{\alpha_c z F}{RT} \eta )
+            \Big)
+
+    where:
+
+        .. math::
+            \eta = V_s - V_{\ell} - V_{oc}
+
+    where ``V_s`` is the solid voltage, ``V_l`` is the electrolyte voltage,
+    and ``V_oc`` is the open-circuit voltage.
+
+    The slope and intercept provide a linearized source term equation
+    about the current value of X as follow:
+
+        .. math::
+            rate = S_{1} X + S_{2}
+
+    Parameters
+    ----------
+    target : GenericPhysics
+        The Physics object where the result will be applied.
+    X : str
+        The dictionary key of the quantity of interest (i.e. main variable
+        to be solved; in this case, concentration).
+    z : float
+        Number of electrons transferred in the redox reaction.
+    j0 : float
+        Exchange current density [A/m^2].
+    c_ref : float
+        Reference concentration [mol/m^3].
+    alpha_anode : float
+        Anodic transfer coefficient.
+    alpha_cathode : float
+        Cathodic transfer coefficient.
+    reaction_order : float
+        Reaction order, i.e. the exponent of the concentration term
+    solid_voltage : str
+        The dictionary key of solid phase voltages [V].
+    electrolyte_voltage : str
+        The dictionary key of electrolyte phase voltages [V].
+    open_circuit_voltage : str
+        The dictionary key of open-circuit voltage values [V].
+    reaction_area : str
+        The dictionary key of reaction area values [m^2].
+    temperature : str
+        The dictionary key of temperature values [K].
+
+    Returns
+    -------
+    dict
+        Dictionary containing the following key/value pairs:
+
+        - rate : The value of the source term function at the given X.
+        - S1 : The slope of the source term function at the given X.
+        - S2 : The intercept of the source term function at the given X.
+
+    Notes
+    -----
+    The difference between butler_volmer_conc and butler_volmer_voltage is
+    that the former is linearized with respect to the electrolyte
+    concentration whereas the latter is linearized with respect to the
+    electrolyte voltage.
+
+    Consequently, while the S1 and S2 produced by these model shouldn't be
+    identical, they should both produce the same **rate** with the only
+    difference that the rate generated by butler_volmer_conc has the units
+    [mol/s] whereas that generated by butler_volmer_voltage has the units
+    [C/s]. Therefore, the two rates will differ by z * F, where z is the
+    number of electrons transferred and F is the Faraday's constant.
+
+    """
+    network = target.project.network
+    pores = network.map_pores(pores=target.Ps, origin=target)
+
+    # Fetch model variables
+    X = target[X]
+    T = target[temperature]
+    Vs = target[solid_voltage]
+    Ve = target[electrolyte_voltage]
+    Voc = target[open_circuit_voltage]
+    A_rxn = network[reaction_area][pores]
+    F = _sp.constants.physical_constants["Faraday constant"][0]
+    R = _sp.constants.R
+    nu = reaction_order
+
+    # Linearize with respect to X (electrolyte concentration)
+    eta = Vs - Ve - Voc
+    cte = j0 * A_rxn / (z * F)
+    m1 = alpha_anode * z * F / (R * T)
+    m2 = alpha_cathode * z * F / (R * T)
+    fV = _np.exp(m1 * eta) - _np.exp(-m2 * eta)
+    fC = (X / c_ref)**nu
+    r = cte * fC * fV
+    drdC = cte * (X / c_ref)**(nu - 1) * (1 / c_ref) * fV
+    S1 = drdC
+    S2 = r - drdC * X
+
+    values = {"pore.S1": S1, "pore.S2": S2, "pore.rate": r}
+    return values
+
+
+def butler_volmer_voltage(
+    target, X, z, j0, c_ref, alpha_anode, alpha_cathode,
+    reaction_order=1,
+    temperature="pore.temperature",
+    reaction_area="pore.reaction_area",
+    solid_voltage="pore.solid_voltage",
+    open_circuit_voltage="pore.open_circuit_voltage",
+    electrolyte_concentration="pore.electrolyte_concentration",
+):
+    r"""
+    Calculates the rate, slope and intercept of the Butler-Volmer kinetic model
+    based on **voltage** to be used in electron conduction algorithms.
+
+        .. math::
+            r = j_0 A_{rxn} (\frac{ c }{ c_{ref} }) ^ {\nu}
+            \Big(
+                \exp(  \frac{\alpha_a z F}{RT} \eta )
+              - \exp( -\frac{\alpha_c z F}{RT} \eta )
+            \Big)
+
+    where:
+
+        .. math::
+            \eta = V_s - X - V_{oc}
+
+    where ``V_s`` is the solid voltage, ``X`` is the electrolyte voltage,
+    and ``V_oc`` is the open-circuit voltage.
+
+    The slope and intercept provide a linearized source term equation
+    about the current value of X as follow:
+
+        .. math::
+            rate = S_{1} X + S_{2}
+
+    Parameters
+    ----------
+    target : GenericPhysics
+        The Physics object where the result will be applied.
+    X : str
+        The dictionary key of the quantity of interest (i.e. main variable
+        to be solved; in this case, electrolyte voltage).
+    z : float
+        Number of electrons transferred in the redox reaction.
+    j0 : float
+        Exchange current density [A/m^2].
+    c_ref : float
+        Reference concentration [mol/m^3].
+    alpha_anode : float
+        Anodic transfer coefficient.
+    alpha_cathode : float
+        Cathodic transfer coefficient.
+    electrolyte_concentration : str
+        The dictionary key of the electrolyte concentrations [mol/m^3].
+    reaction_order : float
+        Reaction order, i.e. the exponent of the concentration term
+    solid_voltage : str
+        The dictionary key of solid phase voltages [V].
+    electrolyte_voltage : str
+        The dictionary key of electrolyte phase voltages [V].
+    open_circuit_voltage : str
+        The dictionary key of open-circuit voltage values [V].
+    reaction_area : str
+        The dictionary key of reaction area values [m^2].
+    temperature : str
+        The dictionary key of temperature values [K].
+
+    Returns
+    -------
+    dict
+        Dictionary containing the following key/value pairs:
+
+        - rate : The value of the source term function at the given X.
+        - S1 : The slope of the source term function at the given X.
+        - S2 : The intercept of the source term function at the given X.
+
+    Notes
+    -----
+    The difference between butler_volmer_conc and butler_volmer_voltage is
+    that the former is linearized with respect to the electrolyte
+    concentration whereas the latter is linearized with respect to the
+    electrolyte voltage.
+
+    Consequently, while the S1 and S2 produced by these model shouldn't be
+    identical, they should both produce the same **rate** with the only
+    difference that the rate generated by butler_volmer_conc has the units
+    [mol/s] whereas that generated by butler_volmer_voltage has the units
+    [C/s]. Therefore, the two rates will differ by z * F, where z is the
+    number of electrons transferred and F is the Faraday's constant.
+
+    """
+    network = target.project.network
+    pores = network.map_pores(pores=target.Ps, origin=target)
+
+    # Fetch model variables
+    A_rxn = network[reaction_area][pores]
+    Vs = target[solid_voltage]
+    Voc = target[open_circuit_voltage]
+    c = target[electrolyte_concentration]
+    T = target[temperature]
+    X = target[X]
+    F = _sp.constants.physical_constants["Faraday constant"][0]
+    R = _sp.constants.R
+    nu = reaction_order
+
+    # Linearize with respect to X (electrolyte voltage)
+    eta = Vs - X - Voc
+    cte = j0 * A_rxn
+    m1 = alpha_anode * z * F / (R * T)
+    m2 = alpha_cathode * z * F / (R * T)
+    fV = _np.exp(m1 * eta) - _np.exp(-m2 * eta)
+    dfVdV = -(m1 * _np.exp(m1 * eta) + m2 * _np.exp(-m2 * eta))
+    fC = (c / c_ref)**nu
+    r = cte * fC * fV
+    drdV = cte * fC * dfVdV
+    S1 = drdV
+    S2 = r - drdV * X
+
+    values = {"pore.S1": S1, "pore.S2": S2, "pore.rate": r}
+    return values


Test patch:
diff --git a/tests/unit/models/physics/GenericSourceTermTest.py b/tests/unit/models/physics/GenericSourceTermTest.py
index 39e59b501..111fb7c0b 100644
--- a/tests/unit/models/physics/GenericSourceTermTest.py
+++ b/tests/unit/models/physics/GenericSourceTermTest.py
@@ -37,7 +37,7 @@ class GenericSourceTermTest:
         self.phase['pore.mole_fraction'] = 0.1
         for source in sources:
             self.phys.add_model(propname="pore.source_term",
-                                model=getattr(pm.generic_source_term, source),
+                                model=getattr(pm.source_terms, source),
                                 X="pore.mole_fraction")
             assert self.phys["pore.source_term.rate"].mean() == 0
             assert self.phys["pore.source_term.S1"].mean() == 0
@@ -47,13 +47,13 @@ class GenericSourceTermTest:
         self.phys['pore.item1'] = 0.5e-11
         self.phys['pore.item2'] = 1.5e-12
         self.phys.add_model(propname='pore.source1',
-                            model=pm.generic_source_term.linear,
+                            model=pm.source_terms.linear,
                             A1='pore.item1',
                             A2='pore.item2',
                             X='pore.mole_fraction',
                             regen_mode='normal')
         self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.linear_sym,
+                            model=pm.source_terms.linear_sym,
                             A1='pore.item1',
                             A2='pore.item2',
                             X='pore.mole_fraction',
@@ -81,14 +81,14 @@ class GenericSourceTermTest:
         self.phys['pore.item2'] = 2.5
         self.phys['pore.item3'] = -1.4e-11
         self.phys.add_model(propname='pore.source1',
-                            model=pm.generic_source_term.power_law,
+                            model=pm.source_terms.power_law,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
                             X='pore.mole_fraction',
                             regen_mode='normal')
         self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.power_law_sym,
+                            model=pm.source_terms.power_law_sym,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -120,7 +120,7 @@ class GenericSourceTermTest:
         self.phys['pore.item5'] = -0.34
         self.phys['pore.item6'] = 2e-14
         self.phys.add_model(propname='pore.source1',
-                            model=pm.generic_source_term.exponential,
+                            model=pm.source_terms.exponential,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -130,7 +130,7 @@ class GenericSourceTermTest:
                             X='pore.mole_fraction',
                             regen_mode='normal')
         self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.exponential_sym,
+                            model=pm.source_terms.exponential_sym,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -164,7 +164,7 @@ class GenericSourceTermTest:
         self.phys['pore.item4'] = -0.34
         self.phys['pore.item5'] = 2e-14
         self.phys.add_model(propname='pore.source1',
-                            model=pm.generic_source_term.natural_exponential,
+                            model=pm.source_terms.natural_exponential,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -173,7 +173,7 @@ class GenericSourceTermTest:
                             X='pore.mole_fraction',
                             regen_mode='normal')
         self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.natural_exponential_sym,
+                            model=pm.source_terms.natural_exponential_sym,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -207,7 +207,7 @@ class GenericSourceTermTest:
         self.phys['pore.item5'] = 0.133
         self.phys['pore.item6'] = -5.1e-13
         self.phys.add_model(propname='pore.source1',
-                            model=pm.generic_source_term.logarithm,
+                            model=pm.source_terms.logarithm,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -217,7 +217,7 @@ class GenericSourceTermTest:
                             X='pore.mole_fraction',
                             regen_mode='normal')
         self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.logarithm_sym,
+                            model=pm.source_terms.logarithm_sym,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -252,7 +252,7 @@ class GenericSourceTermTest:
         self.phys['pore.item4'] = 0.133
         self.phys['pore.item5'] = -5.1e-14
         self.phys.add_model(propname='pore.source1',
-                            model=pm.generic_source_term.natural_logarithm,
+                            model=pm.source_terms.natural_logarithm,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -261,7 +261,7 @@ class GenericSourceTermTest:
                             X='pore.mole_fraction',
                             regen_mode='on_demand')
         self.phys.add_model(propname='pore.source2',
-                            model=pm.generic_source_term.natural_logarithm_sym,
+                            model=pm.source_terms.natural_logarithm_sym,
                             A1='pore.item1',
                             A2='pore.item2',
                             A3='pore.item3',
@@ -298,7 +298,7 @@ class GenericSourceTermTest:
         phys['pore.item4'] = 0.133
         phys['pore.item5'] = -5.1e-14
         phys.add_model(propname='pore.source1',
-                       model=pm.generic_source_term.natural_logarithm,
+                       model=pm.source_terms.natural_logarithm,
                        A1='pore.item1',
                        A2='pore.item2',
                        A3='pore.item3',
@@ -313,13 +313,51 @@ class GenericSourceTermTest:
                                          ('e', 'pore.item5'),
                                          ('x', 'pore.mole_fraction')])
         phys.add_model(propname='pore.general',
-                       model=op.models.physics.generic_source_term.general_symbolic,
+                       model=op.models.physics.source_terms.general_symbolic,
                        eqn=y, arg_map=arg_map,
                        regen_mode='normal')
         assert np.allclose(phys['pore.source1.rate'], phys['pore.general.rate'])
         assert np.allclose(phys['pore.source1.S1'], phys['pore.general.S1'])
         assert np.allclose(phys['pore.source1.S2'], phys['pore.general.S2'])
 
+    def test_butler_volmer_kinetics(self):
+        np.random.seed(10)
+        self.net["pore.reaction_area"] = np.random.rand(self.net.Np)
+        self.phys['pore.electrolyte_voltage'] = np.random.rand(self.net.Np) * 0.2
+        self.phys['pore.solid_voltage'] = 1.1
+        self.phys['pore.open_circuit_voltage'] = 1.2
+        self.phase['pore.electrolyte_concentration'] = np.random.rand(self.net.Np)
+        BV_params = {
+            "z": 4,
+            "j0": 1e-3,
+            "c_ref": 1000,
+            "alpha_anode": 0.4,
+            "alpha_cathode": 0.6
+        }
+        self.phys.add_model(propname='pore.rxn_BV_c',
+                            model=pm.source_terms.butler_volmer_conc,
+                            X="pore.electrolyte_concentration", **BV_params)
+        self.phys.add_model(propname='pore.rxn_BV_v',
+                            model=pm.source_terms.butler_volmer_voltage,
+                            X="pore.electrolyte_voltage", **BV_params)
+        # Check Butler-Volmer model (concentration)
+        S1_BV_c = self.phys["pore.rxn_BV_c.S1"]
+        S2_BV_c = self.phys["pore.rxn_BV_c.S2"]
+        rate_BV_c = self.phys["pore.rxn_BV_c.rate"]
+        eps = np.finfo(float).eps * np.abs(rate_BV_c).max()
+        assert_allclose(S2_BV_c.mean(), 0, atol=eps)
+        assert_allclose(S1_BV_c.mean(), -0.06977055)
+        assert_allclose(rate_BV_c.mean(), -0.03541646)
+        # Check Butler-Volmer model (voltage)
+        S1_BV_v = self.phys["pore.rxn_BV_v.S1"]
+        S2_BV_v = self.phys["pore.rxn_BV_v.S2"]
+        rate_BV_v = self.phys["pore.rxn_BV_v.rate"]
+        assert_allclose(S2_BV_v.mean(), 226884.32)
+        assert_allclose(S1_BV_v.mean(), -1277463.52)
+        assert_allclose(rate_BV_v.mean(), -13668.675)
+        # The two Butler-Volmer models must only differ by z*F (unit conversion)
+        assert_allclose(rate_BV_v, rate_BV_c * BV_params["z"] * 96485.33212)
+
 
 if __name__ == '__main__':
 
diff --git a/tests/unit/phases/mixtures/GenericMixtureTest.py b/tests/unit/phases/mixtures/GenericMixtureTest.py
index 7e7edd1fa..8d4dd2ca3 100644
--- a/tests/unit/phases/mixtures/GenericMixtureTest.py
+++ b/tests/unit/phases/mixtures/GenericMixtureTest.py
@@ -1,6 +1,4 @@
-import pytest
 import numpy as np
-import scipy as sp
 import openpnm as op
 from openpnm.phases import mixtures