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 0x7f2e8c7522a0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Properties Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 pore.coords 9 / 9
3 throat.conns 12 / 12
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Labels Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 pore.xmin 3
3 pore.xmax 3
4 pore.ymin 3
5 pore.ymax 3
6 pore.surface 8
7 throat.surface 8
8 pore.left 3
9 pore.right 3
10 pore.front 3
11 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:
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 theFickianDiffusionalgorithm will call this model to retrieve the reaction rate at the given conditions. However, we still have to tellfdabout this reaction.The model uses the values of
'pore.concentration'to compute \(r_A\), so we specify to use that array forX.Since
'pore.concentration'does not exist yet this model will throw an error if run, so we set theregen_modeto'deferred'.A1,A2andA3can be either scalars or dictionary keys to numpy arrays on thephobject. This is possible since all OpenPNM dict object return any non-string keys directly back, sofd[1.0]returns1.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 0x7f2e8c753010>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# 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 0x7f2e8c753010>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# 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 0x7f2e8c752b60>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# 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:
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:
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:
We could then group all the linear terms as follows:
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:
This can be rearranged to solve for \(S\):
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:
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.40548397 -0.37952706 0. -0.02595691 0. 0.
0. 0. 0. ]
[-0.37952706 1.34886612 -0.55560702 0. -0.41373204 0.
0. 0. 0. ]
[ 0. -0.55560702 0.73380616 0. 0. -0.17819914
0. 0. 0. ]
[-0.02595691 0. 0. 1.45612095 -0.60783281 0.
-0.82233124 0. 0. ]
[ 0. -0.41373204 0. -0.60783281 2.10344055 -0.66855872
0. -0.41331698 0. ]
[ 0. 0. -0.17819914 0. -0.66855872 1.38523868
0. 0. -0.53848081]
[ 0. 0. 0. -0.82233124 0. 0.
1.55205959 -0.72972835 0. ]
[ 0. 0. 0. 0. -0.41331698 0.
-0.72972835 1.93487003 -0.7918247 ]
[ 0. 0. 0. 0. 0. -0.53848081
0. -0.7918247 1.33030551]]