Tutorial 3 of 3: Advanced Topics and Usage

Learning Outcomes

• Use different methods to add boundary pores to a network
• Manipulate network topology by adding and removing pores and throats
• Explore the ModelsDict design, including copying models between objects, and changing model parameters
• Write a custom pore-scale model and a custom Phase
• Access and manipulate objects associated with the network
• Combine multiple algorithms to predict relative permeability

Build and Manipulate Network Topology

For the present tutorial, we'll keep the topology simple to help keep the focus on other aspects of OpenPNM.

Adding Boundary Pores

When performing transport simulations it is often useful to have 'boundary' pores attached to the surface(s) of the network where boundary conditions can be applied. When using the Cubic class, two methods are available for doing this: add_boundaries, which is specific for the Cubic class, and add_boundary_pores, which is a generic method that can also be used on other network types and which is inherited from GenericNetwork. The first method automatically adds boundaries to ALL six faces of the network and offsets them from the network by 1/2 of the value provided as the network spacing. The second method provides total control over which boundary pores are created and where they are positioned, but requires the user to specify to which pores the boundary pores should be attached to. Let's explore these two options:

Adding and Removing Pores and Throats

OpenPNM uses a list-based data storage scheme <topology>_ for all properties, including topological connections. One of the benefits of this approach is that adding and removing pores and throats from the network is essentially as simple as adding or removing rows from the data arrays. The one exception to this 'simplicity' is that the 'throat.conns' array must be treated carefully when trimming pores, so OpenPNM provides the extend and trim functions for added and removing, respectively. To demonstrate, let's reduce the coordination number of the network to create a more random structure:

When the trim function is called, it automatically checks the health of the network afterwards, so logger messages might appear on the command line if problems were found such as isolated clusters of pores or pores with no throats. This health check is performed by calling the Network's' check_network_health method which returns a HealthDict containing the results of the checks:

The HealthDict contains several lists including things like duplicate throats and isolated pores, but also a suggestion of which pores to trim to return the network to a healthy state. Also, the HealthDict has a health attribute that is False is any checks fail.

Define Geometry Objects

The boundary pores we've added to the network should be treated a little bit differently. Specifically, they should have no volume or length (as they are not physically representative of real pores). To do this, we create two separate Geometry objects, one for internal pores and one for the boundaries:

The Stick_and_Ball class is preloaded with the pore-scale models to calculate all the necessary size information (pore diameter, pore.volume, throat lengths, throat.diameter, etc). The GenericGeometry class used for the boundary pores is empty and requires work. Moreover, as throat properties like length and diameter are used in OpenPNM's algorithms (e.g., drainage) so a small value (0.001 nm) has been set to all of them. Furthermore, the volume of throats as cylinderical objects has been calculated from their length and diameter.

These models are required for the Hagan-Poiseuille model. Most of them are straight-forward geometry calculations, except for the model used for 'throat.diameter'. In this case the model looks into the neighbor pores, retrieves the two 'pore.diameter' and uses the 'max' value. Because we set the boundary pores to have 0 diameter, this will naturally find result in the throat being assigned the diameter of the internal pore.

Define Multiple Phase Objects

In order to simulate relative permeability of air through a partially water-filled network, we need to create each Phase object. OpenPNM includes pre-defined classes for each of these common fluids:

Aside: Creating a Custom Phase Class

In many cases you will want to create your own fluid, such as an oil or brine, which may be commonly used in your research. OpenPNM cannot predict all the possible scenarios, but luckily it is easy to create a custom Phase class as follows:

• Creating a Phase class basically involves placing a series of self.add_model commands within the __init__ section of the class definition. This means that when the class is instantiated, all the models are added to itself (i.e. self).
• **kwargs is a Python trick that captures all arguments in a dict called kwargs and passes them to another function that may need them. In this case they are passed to the __init__ method of Oil's parent by the super function. Specifically, things like name and network are expected.
• The above code block also stores the molecular weight of the oil as a constant value
• Adding models and constant values in this way could just as easily be done in a run script, but the advantage of defining a class is that it can be saved in a file (i.e. 'my_custom_phases') and reused in any project.

Define Physics Objects for Each Geometry and Each Phase

In the previous tutorial <intermediate_usage>_ we created two Physics object, one for each of the two Geometry objects used to handle the stratified layers. In this tutorial, the internal pores and the boundary pores each have their own Geometry, but there are two Phases, which also each require a unique Physics:

• To reiterate, one Physics object is required for each Geometry AND each Phase, so the number can grow to become annoying very quickly Some useful tips for easing this situation are given below.

Create a Custom Pore-Scale Physics Model

Perhaps the most distinguishing feature between pore-network modeling papers is the pore-scale physics models employed. Accordingly, OpenPNM was designed to allow for easy customization in this regard, so that you can create your own models to augment or replace the ones included in the OpenPNM models libraries. For demonstration, let's implement the capillary pressure model proposed by Mason and Morrow in 1994. They studied the entry pressure of non-wetting fluid into a throat formed by spheres, and found that the converging-diverging geometry increased the capillary pressure required to penetrate the throat. As a simple approximation they proposed $P_c = -2 \sigma \cdot cos(2/3 \theta) / R_t$

Pore-scale models are written as basic function definitions:

Let's examine the components of above code:

• The function receives network, phase objects as arguments. Each of these provide access to the properties necessary for the calculation: 'pore.diameter' values are retrieved via the network, and the thermophysical properties are retrieved directly from the phase.
• The f value is a scale factor that is applied to the contact angle. Mason and Morrow suggested a value of 2/3 as a decent fit to the data, but we'll make this an adjustable parameter with 2/3 as the default.
• Note the pore.diameter is actually a Geometry property, but it is retrieved via the network using the data exchange rules outlined in the second tutorial, and explained fully in :ref:data_storage.
• All of the calculations are done for every throat in the network, but this pore-scale model is meant to be assigned to a single Physics object. As such, the last line extracts values from the Pc array for the location of physics and returns just the subset.
• The actual values of the contact angle, surface tension, and throat diameter are NOT sent in as numerical arrays, but rather as dictionary keys to the arrays. There is one very important reason for this: if arrays had been sent, then re-running the model would use the same arrays and hence not use any updated values. By having access to dictionary keys, the model actually looks up the current values in each of the arrays whenever it is run.
• It is good practice to include the dictionary keys as arguments, such as sigma = 'throat.contact_angle'. This way the user can control where the contact angle could be stored on the Phase object.

Copy Models Between Physics Objects

As mentioned above, the need to specify a separate Physics object for each Geometry and Phase can become tedious. It is possible to copy the pore-scale models assigned to one object onto another object. First, let's assign the models we need to phys_water_internal:

Now make a copy of the models on phys_water_internal and apply it all the other water Physics objects:

The only 'gotcha' with this approach is that each of the Physics objects must be regenerated in order to place numerical values for all the properties into the data arrays:

Adjust Pore-Scale Model Parameters

The pore-scale models are stored in a ModelsDict object that is itself stored under the models attribute of each object. This arrangement is somewhat convoluted, but it enables integrated storage of models on the object's wo which they apply. The models on an object can be inspected with print(phys_water_internal), which shows a list of all the pore-scale properties that are computed by a model, and some information about the model's regeneration mode.

Each model in the ModelsDict can be individually inspected by accessing it using the dictionary key corresponding to pore-property that it calculates, i.e. print(phys_water_internal)['throat.capillary_pressure']). This shows a list of all the parameters associated with that model. It is possible to edit these parameters directly:

More details about the ModelsDict and ModelWrapper classes can be found in :ref:models.

Perform Multiphase Transport Simulations

Use the Built-In Drainage Algorithm to Generate an Invading Phase Configuration

• The inlet pores were set to both 'top' and 'bottom' using the pn.pores method. The algorithm applies to the entire network so the mapping of network pores to the algorithm pores is 1-to-1.
• The run method automatically generates a list of 25 capillary pressure points to test, but you can also specify more pores, or which specific points to tests. See the methods documentation for the details.
• Once the algorithm has been run, the resulting capillary pressure curve can be viewed with plot_drainage_curve. If you'd prefer a table of data for plotting in your software of choice you can use get_drainage_data which prints a table in the console.

Set Pores and Throats to Invaded

After running, the mip object possesses an array containing the pressure at which each pore and throat was invaded, stored as 'pore.inv_Pc' and 'throat.inv_Pc'. These arrays can be used to obtain a list of which pores and throats are invaded by water, using Boolean logic:

The resulting Boolean masks can be used to manually adjust the hydraulic conductivity of pores and throats based on their phase occupancy. The following lines set the water filled throats to near-zero conductivity for air flow:

• The logic of these statements implicitly assumes that transport between two pores is only blocked if the throat is filled with the other phase, meaning that both pores could be filled and transport is still permitted. Another option would be to set the transport to near-zero if either or both of the pores are filled as well.
• The above approach can get complicated if there are several Geometry objects, and it is also a bit laborious. There is a pore-scale model for this under Physics.models.multiphase called conduit_conductance. The term conduit refers to the path between two pores that includes 1/2 of each pores plus the connecting throat.

Calculate Relative Permeability of Each Phase

We are now ready to calculate the relative permeability of the domain under partially flooded conditions. Instantiate an StokesFlow object:

The relative permeability is the ratio of the water flow through the partially water saturated media versus through fully water saturated media; hence we need to find the absolute permeability of water. This can be accomplished by regenerating the phys_water_internal object, which will recalculate the 'throat.hydraulic_conductance' values and overwrite our manually entered near-zero values from the inv simulation using phys_water_internal.models.regenerate(). We can then re-use the water_flow algorithm:

And finally, the relative permeability can be found from:

• The ratio of the flow rates gives the normalized relative permeability since all the domain size, viscosity and pressure differential terms cancel each other.
• To generate a full relative permeability curve the above logic would be placed inside a for loop, with each loop increasing the pressure threshold used to obtain the list of invaded throats (Ti).
• The saturation at each capillary pressure can be found be summing the pore and throat volume of all the invaded pores and throats using Vp = geom['pore.volume'][Pi] and Vt = geom['throat.volume'][Ti].