

In [1]:

    
import sys
sys.path.append('C:\\Users\\amin\\Documents\\Repos\\OpenPNM')

Relative Diffusivity

Generating the Network, adding Geometry and creating Phases

This example shows you how to calculate a transport property relative to the saturation of the domain by a particular phase. In this case the property is the diffusivity of air relative to the saturation of water. Start by importing OpenPNM and some other useful packages:



In [2]:

    
import openpnm
import numpy as np
import matplotlib.pyplot as plt

Next create a Network object with a cubic topology and lattice spacing of 25 microns and add boundary pores



In [3]:

    
pn = openpnm.network.Cubic(shape=[10, 10, 10], spacing=4.0e-5)
pn.add_boundary_pores()
proj = pn.project

Next create a Geometry to manage the pore and throat size information. A Geometry can span over a part of the Network only, so we need to specify to which pores and throats this Geometry object should apply. For this example, we want one Geometry to apply to the internal pores and a different one to the boundary pores:



In [4]:

    
Ps = pn.pores('*boundary', mode='not')
Ts = pn.find_neighbor_throats(pores=Ps, mode='intersection', flatten=True)
geom = openpnm.geometry.StickAndBall(network=pn, pores=Ps, throats=Ts)
Ps = pn.pores('*boundary')
Ts = pn.find_neighbor_throats(pores=Ps, mode='exclusive_or')

The Toray090 Geometry is a predefined class that applies Weibull distributed pore and throat sizes to the internal pores and is named after the Toray carbon paper gas diffusion layer commonly used in fuel cells. The Boundary class is predefined with properties suitable for boundaries such as 0 volume and length.

We must also create the Phase objects, for our purposes the standard air and water phase classes provided are fine:



In [5]:

    
air = openpnm.phases.Air(network=pn, name='air')
water = openpnm.phases.Water(network=pn, name='water')

Define the Pore-Scale Physics

For this simulation the standard physics object can be used as it contains capillary pressure for use in the OrdinaryPercolation algorithm and diffusive conductance for use in the FickianDiffusion algorithm.



In [6]:

    
phys_air = openpnm.physics.Standard(network=pn, phase=air, geometry=geom)
phys_water = openpnm.physics.Standard(network=pn, phase=water, geometry=geom)









    



2018-07-17 12:23:43,359 | WARNING  | root._regen | throat.entry_pressure was not run since the following properties are missing: ['pore.surface_tension', 'throat.diameter']
2018-07-17 12:23:43,377 | WARNING  | root._regen | throat.electrical_conductance was not run since the following properties are missing: ['pore.electrical_conductivity', 'throat.electrical_conductivity', 'throat.equivalent_area', 'throat.conduit_lengths']

Set up and run the Ordinary Percolation Algorithm

In order to simulate a partially saturated material we first run an OrdinaryPercolation Algorithm which sequentially invades the network with the invading_phase based on the capillary pressure of the throats in the network. This allows us to inspect which pores and throats are occupied by which phase at various capillary pressures and this occupancy is used to calculate the multiphase diffusive conductance.



In [11]:

    
bounds = [['front', 'back'], ['left', 'right'], ['top', 'bottom']]
OP_1 = openpnm.algorithms.OrdinaryPercolation(network=pn,
                                              invading_phase=water,
                                              defending_phase=air)
inlets = pn.pores('bottom_boundary')
step = 2
used_inlets = [inlets[x] for x in range(0, len(inlets), step)]
OP_1.run(inlets=used_inlets)
OP_1.return_results()









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-e66c3e2e5311> in <module>()
      2 OP_1 = openpnm.algorithms.OrdinaryPercolation(network=pn,
      3                                               invading_phase=water,
----> 4                                               defending_phase=air)
      5 inlets = pn.pores('bottom_boundary')
      6 step = 2

~\Documents\Repos\OpenPNM\openpnm\algorithms\OrdinaryPercolation.py in __init__(self, settings, **kwargs)
     78 
     79     def __init__(self, settings={}, **kwargs):
---> 80         super().__init__(**kwargs)
     81         self.settings.update(default_settings)
     82         # Use the reset method to initialize all arrays

~\Documents\Repos\OpenPNM\openpnm\algorithms\GenericAlgorithm.py in __init__(self, network, project, settings, **kwargs)
     23         if network is not None:
     24             project = network.project
---> 25         super().__init__(project=project, **kwargs)
     26 
     27         if project.network is not None:

TypeError: __init__() got an unexpected keyword argument 'invading_phase'

Here we have selected half of the boundary pores at the bottom of the domain as inlets for the percolation Algorithm. OrdinaryPercolation has a helpful plotting function which displays the saturation of the invading phase (volume fraction of the pore space) vs. capillary pressure: OP_1.plot_drainage_curve(). The red line is pore saturation, blue is throat saturation and green is total saturation.

Note : OpenPNM's Postprocessing module also has the ability to plot drainage curves and is suitable for use with the InvasionPercolation algorithm too.

Run a Fickian Diffusion Algorithm for each step of the invasion process

We now need to model how the presence of the phases affects the diffusive conductivity of the network. Currently the Physics objects have a property called throat.diffusive_conductance but this model does not account for the occupancy of each phase and assumes that the phase occupies every pore-throat-pore conduit. OpenPNM has a number of multiphase models including a conduit conductance that multiplies the single phase conductance by a factor (default 0.000001) when the phase associated with the physics object is not present. The model has a mode which defaults to 'strict' which applies the conductivity reduction if any one of the connected pores or connecting throat is unoccupied.



In [8]:

    
import openpnm.models.physics as pm
phys_air.add_model(model=pm.multiphase.conduit_conductance,
                   propname='throat.conduit_diffusive_conductance',
                   throat_conductance='throat.diffusive_conductance')
phys_water.add_model(model=pm.multiphase.conduit_conductance,
                     propname='throat.conduit_diffusive_conductance',
                     throat_conductance='throat.diffusive_conductance')









    



2018-07-17 12:23:57,476 | WARNING  | root._regen | throat.conduit_diffusive_conductance was not run since the following properties are missing: ['throat.occupancy', 'pore.occupancy']
2018-07-17 12:23:57,479 | WARNING  | root._regen | throat.conduit_diffusive_conductance was not run since the following properties are missing: ['throat.occupancy', 'pore.occupancy']

Now we create some variables to store our data in for each principle direction (x, y, z). The boundary planes at each side of the domain are used as boundary pores for the Diffusion algorithm.



In [9]:

    
bounds = [['front', 'back'], ['left', 'right'], ['top', 'bottom']]
diff_air = {'0': [], '1': [], '2': []}
diff_water = {'0': [], '1': [], '2': []}
sat= []
tot_vol = np.sum(pn["pore.volume"]) + np.sum(pn["throat.volume"])

Now for each invasion step we cycle through the principle directions and create FickianDiffusion objects for each phase and calculate the effective diffusivity.



In [10]:

    
for Pc in np.unique(OP_1['pore.inv_Pc']):
    OP_1.return_results(Pc=Pc)
    phys_air.regenerate()
    phys_water.regenerate()
    this_sat = 0
    this_sat += np.sum(pn["pore.volume"][water["pore.occupancy"] == 1])
    this_sat += np.sum(pn["throat.volume"][water["throat.occupancy"] == 1])
    sat.append(this_sat)
    for bound_increment in range(len(bounds)):
        BC1_pores = pn.pores(labels=bounds[bound_increment][0]+'_boundary')
        BC2_pores = pn.pores(labels=bounds[bound_increment][1]+'_boundary')
        FD_1 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
        FD_1.set_boundary_conditions(bctype='Dirichlet', bcvalue=0.6,
                                     pores=BC1_pores)
        FD_1.set_boundary_conditions(bctype='Dirichlet', bcvalue=0.2,
                                     pores=BC2_pores)
        FD_1.run(conductance='conduit_diffusive_conductance')
        eff_diff = FD_1.calc_eff_diffusivity()
        diff_air[str(bound_increment)].append(eff_diff)
        FD_2 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=water)
        FD_2.set_boundary_conditions(bctype='Dirichlet', bcvalue=0.6,
                                     pores=BC1_pores)
        FD_2.set_boundary_conditions(bctype='Dirichlet', bcvalue=0.2,
                                     pores=BC2_pores)
        FD_2.run(conductance='conduit_diffusive_conductance')
        eff_diff = FD_2.calc_eff_diffusivity()
        diff_water[str(bound_increment)].append(eff_diff)
        workspace.purge_object(FD_1)
        workspace.purge_object(FD_2)









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-4d6df75a253a> in <module>()
----> 1 for Pc in np.unique(OP_1['pore.inv_Pc']):
      2     OP_1.return_results(Pc=Pc)
      3     phys_air.regenerate()
      4     phys_water.regenerate()
      5     this_sat = 0

NameError: name 'OP_1' is not defined

The return_results method updates the two Phase objects with the occupancy at the given capillary pressure (Pc). The Physics objects are then regenerated to re-calculate the conduit_diffusive_conductance property.

Note : Six Diffusion algorithm objects could have been created outside the loop and then run over and over with the updated conductance values and this would possibly save some computational time but is not as nice to code. Instead the objects are purged and redefined with updated boundary pores inside another for loop.

Plot the Relative Diffusivity Curves for each direction and Phase

Now tidy up the data converting them into Numpy arrays for easy plotting and manipulation and normalize the results by the single phase values:



In [18]:

    
sat=np.asarray(sat)
sat /= tot_vol
rel_diff_air_x    =  np.asarray(diff_air['0'])
rel_diff_air_x   /=  rel_diff_air_x[0]
rel_diff_air_y    =  np.asarray(diff_air['1'])
rel_diff_air_y   /=  rel_diff_air_y[0]
rel_diff_air_z    =  np.asarray(diff_air['2'])
rel_diff_air_z   /=  rel_diff_air_z[0]
rel_diff_water_x  =  np.asarray(diff_water['0'])
rel_diff_water_x /=  rel_diff_water_x[-1]
rel_diff_water_y  =  np.asarray(diff_water['1'])
rel_diff_water_y /=  rel_diff_water_y[-1]
rel_diff_water_z  =  np.asarray(diff_water['2'])
rel_diff_water_z /=  rel_diff_water_z[-1]









    



---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-18-02a093f6b455> in <module>()
      2 sat /= tot_vol
      3 rel_diff_air_x    =  np.asarray(diff_air['0'])
----> 4 rel_diff_air_x   /=  rel_diff_air_x[0]
      5 rel_diff_air_y    =  np.asarray(diff_air['1'])
      6 rel_diff_air_y   /=  rel_diff_air_y[0]

IndexError: index 0 is out of bounds for axis 0 with size 0

Finally plot the results:



In [20]:

    
[lt.ioff()  # Turn off automatic plotting
fig = plt.figure()
ax = fig.gca()
plots=[]
plots.append(plt.plot(sat,rel_diff_air_x,'^-r',label='Dra_x'))
plots.append(plt.plot(sat,rel_diff_air_y,'^-g',label='Dra_y'))
plots.append(plt.plot(sat,rel_diff_air_z,'^-b',label='Dra_z'))
plots.append(plt.plot(sat,rel_diff_water_x,'*-r',label='Drw_x'))
plots.append(plt.plot(sat,rel_diff_water_y,'*-g',label='Drw_y'))
plots.append(plt.plot(sat,rel_diff_water_z,'*-b',label='Drw_z'))
h = plt.xlabel('Liquid Water Saturation')
h = plt.ylabel('Relative Diffusivity')
handles, labels = ax.get_legend_handles_labels()
h = ax.legend(handles, labels,loc=1)
box = ax.get_position()
h = ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])
h = ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))









    



  File "<ipython-input-20-ead7504c308f>", line 2
    fig = plt.figure()
      ^
SyntaxError: invalid syntax

Finally, to draw the plot type fig.canvas.draw_idle() to get the following: