Getting Started#

It is relatively easy to get going with a quick simulation in OpenPNM. In fact the following code block produces a mercury intrusion simulation in just a few lines.

Creating a Cubic Network#

import openpnm as op
import numpy as np
import matplotlib.pyplot as plt
op.visualization.set_mpl_style()

Nx, Ny, Nz = 10, 10, 10
Lc = 1e-4
pn = op.network.Cubic([Nx, Ny, Nz], spacing=Lc)
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()
print(pn)
══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Cubic at 0x7fa8afa1a3e0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.coords                                                   1000 / 1000
  3  throat.conns                                                  2700 / 2700
  4  pore.coordination_number                                      1000 / 1000
  5  pore.max_size                                                 1000 / 1000
  6  throat.spacing                                                2700 / 2700
  7  pore.seed                                                     1000 / 1000
  8  pore.diameter                                                 1000 / 1000
  9  throat.max_size                                               2700 / 2700
 10  throat.diameter                                               2700 / 2700
 11  throat.cross_sectional_area                                   2700 / 2700
 12  throat.hydraulic_size_factors                                 2700 / 2700
 13  throat.diffusive_size_factors                                 2700 / 2700
 14  throat.lens_volume                                            2700 / 2700
 15  throat.length                                                 2700 / 2700
 16  throat.total_volume                                           2700 / 2700
 17  throat.volume                                                 2700 / 2700
 18  pore.volume                                                   1000 / 1000
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.xmin                                                             100
  3  pore.xmax                                                             100
  4  pore.ymin                                                             100
  5  pore.ymax                                                             100
  6  pore.zmin                                                             100
  7  pore.zmax                                                             100
  8  pore.surface                                                          488
  9  throat.surface                                                        972
 10  pore.left                                                             100
 11  pore.right                                                            100
 12  pore.front                                                            100
 13  pore.back                                                             100
 14  pore.bottom                                                           100
 15  pore.top                                                              100
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Defining a Phase#

hg = op.phase.Mercury(network=pn)
hg.add_model(propname='throat.entry_pressure',
             model=op.models.physics.capillary_pressure.washburn)
hg.regenerate_models()
print(hg)
══════════════════════════════════════════════════════════════════════════════
phase_01 : <openpnm.phase.Mercury at 0x7fa85e9f85e0>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.temperature                                              1000 / 1000
  3  pore.pressure                                                 1000 / 1000
  4  throat.contact_angle                                          2700 / 2700
  5  pore.thermal_conductivity                                     1000 / 1000
  6  pore.surface_tension                                          1000 / 1000
  7  pore.viscosity                                                1000 / 1000
  8  pore.density                                                  1000 / 1000
  9  pore.molar_density                                            1000 / 1000
 10  throat.entry_pressure                                         2700 / 2700
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                             1000
  3  throat.all                                                           2700
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Performing a Drainage Simulation#

mip = op.algorithms.Drainage(network=pn, phase=hg)
mip.set_inlet_BC(pores=pn.pores(['left', 'right']))
mip.run(pressures=np.logspace(4, 6))
data = mip.pc_curve()
fig, ax = plt.subplots(figsize=(5.5, 4))
ax.semilogx(data.pc, data.snwp, 'k-o')
ax.set_xlabel('capillary pressure [Pa]')
ax.set_ylabel('mercury saturation');
../_images/f96f84c9a1b924ebd15ac28666ec99d6fb1c4e291d12aba2f7e87e35923451c5.png

Calculating Permeability Coefficient#

As another example, the permeability coefficient can be found as follows:

# Generate phase and physics
water = op.phase.Water(network=pn)
water.add_model(propname='throat.hydraulic_conductance',
                model=op.models.physics.hydraulic_conductance.generic_hydraulic)

# Create algorithm, set boundary conditions and run simulation
sf = op.algorithms.StokesFlow(network=pn, phase=water)
Pin, Pout = (200_000, 101_325)
sf.set_value_BC(pores=pn.pores('left'), values=Pin)
sf.set_value_BC(pores=pn.pores('right'), values=Pout)
sf.run()

The total flow rate into the domain through the boundary pores can be found using sf.rate(pores=pn.pores('xmin')). The permeability coefficient can be found by inserting known values into Darcy’s law as follows:

Q = sf.rate(pores=pn.pores('left'))
A = Ny*Nz*Lc**2
L = Nx*Lc
mu = water['pore.viscosity'].mean()
K = Q*mu*L/(A*(Pin-Pout))
print(K)
[6.42847014e-13]

Adjusting Pore Size Distribution#

It’s also worth explaining how to adjust the pore size distribution of the network, so that the capillary curve and permeability coefficient can be changed to match known values. The geo object controls the geometric properties, and it possess models to calculate values on demand. Let’s change the pore size distribution to a Weibull distribution, but first let’s store the existing values in a dummy variable so we can compare later.

import openpnm.models.geometry as gmods

pn['pore.old_diameter'] = pn.pop('pore.diameter')
pn.add_model(propname='pore.diameter',
             model=gmods.pore_size.weibull,
             shape=0.5, loc=0, scale=1e-5)

Now you can use matplotlib. to get a quick glance at the histograms of the two distributions.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6, 4.5))
ax.hist(pn['pore.diameter'], edgecolor='k', label='new diameter')
ax.hist(pn['pore.old_diameter'], edgecolor='k', label='old_diameter', bins=20)
ax.set_xlabel('diameter [m]')
ax.set_ylabel('frequency')
ax.legend();
../_images/e7409eb629a15141d4c8b1bff80a0a8fa80839bc28272d1266c94395925e6f3f.png

More complex tasks are explained in the examples page.