Reactive Transport#

This notebook is still a work in progress

Incorporating chemical reactions or source terms in general is an important feature. This notebook will explain how OpenPNM implements this.

import openpnm as op
import numpy as np
op.visualization.set_mpl_style()
pn = op.network.Cubic(shape=[3, 3, 1])
print(pn)
══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Cubic at 0x7f139c737f40>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.coords                                                         9 / 9
  3  throat.conns                                                      12 / 12
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.surface                                                            8
  3  throat.surface                                                          8
  4  pore.left                                                               3
  5  pore.right                                                              3
  6  pore.front                                                              3
  7  pore.back                                                               3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

And let’s also create a Phase object to store the diffusive conductance values:

ph = op.phase.Phase(network=pn)
ph['throat.diffusive_conductance'] = np.random.rand(pn.Nt)

Using the source term functionality in OpenPNM#

First, we’ll illustrate the process of adding a source term to a diffusion simulation, then we’ll look behind the scenes at how it works.

fd = op.algorithms.FickianDiffusion(network=pn, phase=ph)
fd.set_value_BC(pores=pn.pores('left'), values=1.0)

OpenPNM’s models module contains a variety to pre-written source terms. Let’s use a standar reaction, also known as a power-law:

\[ y = A_0 X^{A_1} + A_3 \rightarrow r_A = kC_A^b\]

where \(C_A\) is the concentration of the diffusion species \(A\), \(k\) and \(b\) are the kinetic constant (actually a strong function of temperature) and \(b\) is the reaction order.

f = op.models.physics.source_terms.power_law

This model is added to the ph object, since it is a function of the phase properties, specifically concentration and temperature if the reaction is not isothermal.

ph.add_model(propname='pore.reaction', 
             model=f,
             X='pore.concentration',
             A1=1e-5,
             A2=1,
             A3=0.0,
             regen_mode='deferred')

Let’s unpack the above code block:

  • The model is being stored under 'pore.reaction', which means that the FickianDiffusion algorithm will call this model to retrieve the reaction rate at the given conditions. However, we still have to tell fd about this reaction.

  • The model uses the values of 'pore.concentration' to compute \(r_A\), so we specify to use that array for X.

  • Since 'pore.concentration' does not exist yet this model will throw an error if run, so we set the regen_mode to 'deferred'.

  • A1, A2 and A3 can be either scalars or dictionary keys to numpy arrays on the ph object. This is possible since all OpenPNM dict object return any non-string keys directly back, so fd[1.0] returns 1.0.

Now we add this model to our algorithm, to tell the algorithm about it and also in which pores it applies:

fd.set_source(propname='pore.reaction', pores=pn.pores('right'))

The act of adding a source term creates a new array on the algorithm as can be seen below:

print(fd)
══════════════════════════════════════════════════════════════════════════════
fick_02 : <openpnm.algorithms.FickianDiffusion at 0x7f139c2af4a0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.bc.rate                                                        0 / 9
  3  pore.bc.value                                                       3 / 9
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                                9
  3  throat.all                                                             12
  4  pore.source.reaction                                                    3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

The 'pore.source.reaction' array contains a boolean array where True values indicate the locations where the source term should be applied. Had we added a second reaction, called 'pore.another_reaction', the above printout would have an entry for 'pore.source.another_reaction'. This is designed to make it easy to get a list of all source terms and their location using fd['pore.source'] which would return a dict with 'reaction' and 'another_reaction' as keys, and two boolean arrays as values indicating where each reaction was applied.

So now the algorithm is ready to be run:

fd.run()

Let’s look at which arrays are created on both fd and ph:

print(fd)
══════════════════════════════════════════════════════════════════════════════
fick_02 : <openpnm.algorithms.FickianDiffusion at 0x7f139c2af4a0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.bc.rate                                                        0 / 9
  3  pore.bc.value                                                       3 / 9
  4  pore.concentration                                                  9 / 9
  5  pore.initial_guess                                                  9 / 9
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                                9
  3  throat.all                                                             12
  4  pore.source.reaction                                                    3
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Above we can see that 'pore.concentration' was created, which contains the solution of the transport problem. 'pore.initial_guess' was also created since the algorithm was aware that a source term was added, so iteration of the problem would be required hence an initial guess was needed.

print(ph)
══════════════════════════════════════════════════════════════════════════════
phase_01 : <openpnm.phase.Phase at 0x7f139c2af5e0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.temperature                                                    9 / 9
  3  pore.pressure                                                       9 / 9
  4  throat.diffusive_conductance                                      12 / 12
  5  pore.concentration                                                  9 / 9
  6  pore.reaction.S1                                                    9 / 9
  7  pore.reaction.S2                                                    9 / 9
  8  pore.reaction.rate                                                  9 / 9
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                                9
  3  throat.all                                                             12
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

We can see that 'pore.concentration' was also added to this object. This was necessary since the source-term models needed access to this array in order to compute the reaction rate. So the fd algorithm writes this array to fd and then regenerates the source term model(s). Recall that we set regen_mode='deferred' when assigning this model, so the algorithm is the first to run this model.

We can also see 'pore.reaction.S1', 'S2' and 'rate' which are all values returned by the model. We’ll explain what 'S1' and 'S2' mean in the following section, and 'rate' is self explanatory.

How the source-term functionality works#

OpenPNM uses the method described by Patankar [ref] to linearize the source-terms. As will be shown below, this method computes the slope and intercept of non-linear term at the current value of the quantity being computed. This is the origin of the 'S1' (the slope) and 'S2' (the intercept) arrays on the ph object above.

Determination of 'S1' and 'S2' will be outlined below, but let’s first illustrate how they are used. Consider the mass balance around pore 6, which connects to pore 3 and 7:

\[ g_{3,6} (C_{A,3} - C_{A,6}) + g_{7,6} (C_{A,7} - C_{A,6}) = r_A \]

Meaning the net diffusive flux into pore 6 from pores 3 and 7 is equal to the rate of consumption in pore 6. We can further write the rate value as:

\[ g_{3,6} (C_{A,3} - C_{A,6}) + g_{7,6} (C_{A,7} - C_{A,6}) = k C_{A,6}^2 \]

Now we can see there will be a problem solving this equation since a non-linear concentration term has appeared. If we could just write the reaction term as a linear expresssion we would be able to proceed:

\[ g_{3,6} (C_{A,3} - C_{A,6}) + g_{7,6} (C_{A,7} - C_{A,6}) = S_1 C_{A,6} + S_2 \]

We could then group all the linear terms as follows:

\[ g_{3,6} C_{A,3} + (-g_{3,6} - g_{7,6} - S_1) C_{A,6} + g_{7,6} C_{A,7} = S_2 \]

The above equation (or rather the system of equations to which is belongs) is then solved using standard numerical solvers to find all the unknown values of \(C_A\). Now we just need to find 'S1' and 'S2'.

The process starts with the standard definition of Newton’s method:

\[ S'=\frac{(S-S^*)}{(x_i-x_i^* )} \]

This can be rearranged to solve for \(S\):

\[ S=S^*+(S')^*⋅(x_i-x_i^* ) \]

Where \(S'^*\) is is the derivative of the function at \(x_i^*\). Since \(x_i^*\) is known (as either an initial guess or a updated value from each iteration), we can treat it and \(S'^*\) as constants and lump them together as follows:

\[ S= (S' )^* x_i+(S^*-(S' )^* x_i^* )=S_1 x_i+S_2 \]

So both \(S\) values are computed from known values, and the result is a linear equation that is a function of the unknown \(x_i\), which can be entered into the system of equations as mentioned in the previous section.

Now that we have seen how 'S1' and 'S2' are defined and computed, we can revisit how this is implemented in OpenPNM. We the run method is called on the fd, several hidden methods are called, these include:

  • _build_A

  • _build_b

  • _apply_bcs

  • _apply_sources

Let see how the coefficient and RHS matrix evolve as these are each called.

fd = op.algorithms.FickianDiffusion(network=pn, phase=ph)
fd.set_value_BC(pores=pn.pores('left'), values=1.0)
fd.set_source(propname='pore.reaction', pores=pn.pores('right'))
print(fd.A.todense())
[[ 0.90695301 -0.78004082  0.         -0.12691219  0.          0.
   0.          0.          0.        ]
 [-0.78004082  2.01607438 -0.82785476  0.         -0.40817879  0.
   0.          0.          0.        ]
 [ 0.         -0.82785476  1.38596348  0.          0.         -0.55810872
   0.          0.          0.        ]
 [-0.12691219  0.          0.          1.48171159 -0.55730265  0.
  -0.79749675  0.          0.        ]
 [ 0.         -0.40817879  0.         -0.55730265  2.03183203 -0.14252252
   0.         -0.92382807  0.        ]
 [ 0.          0.         -0.55810872  0.         -0.14252252  1.35979693
   0.          0.         -0.65916569]
 [ 0.          0.          0.         -0.79749675  0.          0.
   1.45602018 -0.65852344  0.        ]
 [ 0.          0.          0.          0.         -0.92382807  0.
  -0.65852344  1.85753327 -0.27518177]
 [ 0.          0.          0.          0.          0.         -0.65916569
   0.         -0.27518177  0.93434746]]