Basic Usage#
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.
Problem setup#
import openpnm as op
# Define geometrical parameters
Lc = 1e-4
Nx, Ny, Nz = (10, 10, 10)
# Generate network, geometry, phase, and physics
pn = op.network.Cubic(shape=[Nx, Ny, Nz], spacing=Lc)
geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)
Hg = op.phases.Mercury(network=pn)
phys = op.physics.Standard(network=pn, phase=Hg, geometry=geo)
# Create algorithm and run simulation
mip = op.algorithms.Porosimetry(network=pn, phase=Hg)
mip.set_inlets(pores=pn.pores(['left', 'right', 'top', 'bottom']))
mip.run()
Basic visualization#
The results can be visualized with mip.plot_intrusion_curve()
giving
something like this:

Calculating permeability#
As another example, the permeability coefficient can be found as follows:
# Generate phase and physics
water = op.phases.Water(network=pn)
phys = op.physics.Standard(network=pn, phase=water, geometry=geo)
# 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('left'))
. 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))
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 op.models.geometry as gmods
geo['pore.old_diameter'] = geo.pop('pore.diameter')
geo.add_model(propname='pore.diameter',
model=gmods.pore_size.weibull,
shape=0.5, loc=0, scale=1e-5)
Now you can run geo.show_hist(['pore.old_diameter', 'pore.diameter'])
to get a quick glance at the histograms of the two distributions.
More complex tasks are explained in the examples page.