Capillary number (flow rate) question

[qingqi-code]

Hi, first of all thank you very much for the developer team effort. I want to use OpenPNM to analyze the Capillary number (Ca or Nc) influence on capillary pressure and relative permeability. Changing flow rate can lead to Capillary number change therefore affect capillary pressure curve and relatively permeability curve (https://link.springer.com/article/10.1007/BF00136349). When I implement flow rate Boundary condition in MixedInvasionPercolation algorithm (I found a function called apply_flow)to calculate capillary pressure and StokesFlow algorithm (function set_rate_BC), the curves have no change.
For the Capillary pressure:

import numpy as np
import pandas as pd
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import openpnm.models.physics as pmods
import matplotlib.pyplot as plt

N = 100
pn = op.network.Cubic(shape=[N, N, 1], spacing=2.5e-5)
geom = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)
water = op.phases.Water(network=pn)
phys = op.physics.Standard(network=pn, phase=water, geometry=geom)
phys['pore.entry_pressure'] = 0.0
fig, ax = plt.subplots(figsize=[5, 5])
ax.hist(phys['throat.entry_pressure'])
plt.show()

def run_mp(trapping=False, residual=None, snap_off=False, plot=True, flowrate=None, phase=None):
    alg = op.algorithms.MixedInvasionPercolation(network=pn)
    if snap_off:
        alg.settings['snap_off'] = 'throat.snap_off'
    alg.setup(phase=phase)
    alg.set_inlets(pores=pn.pores('left'))
    if residual is not None:
        alg.set_residual(pores=residual)
    alg.run()
    if trapping:
        alg.set_outlets(pn.pores('right'))
        alg.apply_trapping()
    inv_points = np.arange(0, 100, 1)
    # returns data as well as plotting
    alg_data = alg.get_intrusion_data(inv_points=inv_points)
    water.update(alg.results(Pc=inv_points.max()))
    if plot:
        fig, ax = plt.subplots(figsize=[5, 5])
        L = np.sqrt(pn.Np).astype(int)
        ax.imshow(alg['pore.invasion_sequence'].reshape([L, L]),
                   cmap=plt.get_cmap('Blues'))
        plt.show()
    if flowrate is not None:
        alg.apply_flow(flowrate=flowrate)
    return alg

alg1 = run_mp(phase=water, flowrate=.1)
alg1_data = alg1.get_intrusion_data()

alg2 = run_mp(phase=water, flowrate=20000)
alg2_data = alg2.get_intrusion_data()

fig, ax = plt.subplots(figsize=[4*1.25, 3*1.25])
ax.plot(alg1_data.S_tot, alg1_data.Pcap);
ax.plot(alg2_data.S_tot, alg2_data.Pcap);

the map and curves are absolutely same:

And the code for Relative Permeability:

import openpnm as op
import numpy as np
import matplotlib.pyplot as plt
import os
from pathlib import Path

INST_DIR = os.path.dirname(os.path.abspath(op.__file__))
ROOT_DIR = Path(INST_DIR).parent
BEREA_DIR = os.path.join(ROOT_DIR, 'OpenPNM-dev\\tests\\fixtures\\ICL-Sandstone(Berea)')
project = op.io.Statoil.import_data(path=BEREA_DIR, prefix='Berea')

pn = project.network
pn['pore.diameter'] = pn['pore.radius'] * 2
pn['throat.diameter'] = pn['throat.radius'] * 2
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
h = pn.check_network_health()

water = op.phases.Water(network=pn, name='water')
oil = op.phases.GenericPhase(network=pn, name='oil')

oil['pore.viscosity'] = 0.001
oil['pore.contact_angle'] = 120
water['pore.contact_angle'] = 60
oil['pore.surface_tension'] = 0.04
water["pore.surface_tension"] = 0.04
water["pore.temperature"] = 350
water.regenerate_models()
oil.regenerate_models()

phys_water = op.physics.GenericPhysics(network=pn, phase=water, geometry=pn)
phys_oil = op.physics.GenericPhysics(network=pn, phase=oil, geometry=pn)
mod = op.models.physics.capillary_pressure.washburn
water.add_model(propname='throat.entry_pressure', model=mod)
oil.add_model(propname='throat.entry_pressure', model=mod)

mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys_water.add_model(propname='throat.hydraulic_conductance', model=mod)
phys_oil.add_model(propname='throat.hydraulic_conductance', model=mod)

IP_withdrawal = op.algorithms.MixedInvasionPercolation(network=pn)
IP_withdrawal.setup(phase=oil)
# Select every 4th pore as inlet to prevent closing it off straight away
IP_withdrawal.set_inlets(pores=pn.pores('inlets')[::4])
IP_withdrawal.run()
IP_withdrawal.set_outlets(pores=pn.pores(['outlets']))
IP_withdrawal.apply_trapping()
injection_data = IP_withdrawal.get_intrusion_data()

Pcdata = {}
Pcdata['capillary_pressure'] = injection_data.Pcap
Pcdata['oil_phase_saturation'] = injection_data.S_tot
Pcdata['water_phase_saturation'] = []
for i in injection_data.S_tot:
    Pcdata['water_phase_saturation'].append(1-i)

FlowAlg_water = op.algorithms.StokesFlow(network=pn, phase=water)
#FlowAlg_water.set_value_BC(pores=pn.pores('inlets'), values=202650)
FlowAlg_water.set_rate_BC(rate=150,pores=pn.pores('inlets'))
FlowAlg_water.set_value_BC(pores=pn.pores('outlets'), values=101325)
FlowAlg_water.run()

FlowAlg_oil = op.algorithms.StokesFlow(network=pn, phase=oil)
#FlowAlg_oil.set_value_BC(pores=pn.pores('inlets'), values=202650)
FlowAlg_oil.set_rate_BC(rate=150,pores=pn.pores('inlets'))
FlowAlg_oil.set_value_BC(pores=pn.pores('outlets'), values=101325)
FlowAlg_oil.run()

Qtot_water = FlowAlg_water.rate(pores=pn.pores('outlets'))
Qtot_oil = FlowAlg_oil.rate(pores=pn.pores('outlets'))

pnVolume = np.sum(pn['pore.volume']) + np.sum(pn['throat.volume'])

Krdata = {}
Krdata['Krw'] = []
Krdata['Krnw'] = []
Krdata['Sw'] = []
for Pc in np.linspace(0, 20000, 21):
    res = IP_withdrawal.results(Pc)
    Po = res['pore.occupancy'] > 0
    To = res['throat.occupancy'] > 0
    oil.update({'pore.occupancy': Po,
                'throat.occupancy': To})
    phys_oil['throat.hydraulic_conductance'][~To] = 1e-20
    FlowAlg_oil.run()
    Q_partial = FlowAlg_oil.rate(pores=pn.pores('outlets'))
    Q_partial[Q_partial < 0] = 0
    Krdata['Krnw'].append(Q_partial[0]/Qtot_oil[0])
    Vt = np.sum(pn['throat.volume'][To])
    Vp = np.sum(pn['pore.volume'][Po])
    Krdata['Sw'].append(1-(Vp+Vt)/pnVolume)
    oil.regenerate_models()
    phys_oil.regenerate_models()
    phys_water['throat.hydraulic_conductance'][To] = 1e-20
    FlowAlg_water.run()
    Q_partial = FlowAlg_water.rate(pores=pn.pores('outlets'))
    Q_partial[Q_partial < 0] = 0
    Krdata['Krw'].append(Q_partial[0]/Qtot_water[0])
    water.regenerate_models()
    phys_water.regenerate_models()

plt.figure()
plt.plot(Krdata['Sw'], Krdata['Krnw'], 'y*-',
         label='Non-wet phase (oil)')
plt.plot(Krdata['Sw'], Krdata['Krw'], 'b*-',
         label='Wet phase (water)')
plt.ylabel('Relative permeability')
plt.xlabel('Saturation \n Porous Volume Fraction occupied by wet phase')
plt.grid(True)
plt.legend()
plt.show()

No matter what flow rate boundary condition I apply (1.5, 150...), the plot is same:

I don't know how to fix with this problem. Will the flow rate boundary condition (or Caplliary number) have influence on Capillary pressure and Relative permeability? Thank you so much!

[jgostick]

Hi @qingqi-code
I'm sorry, but OpenPNM does not have a 'dynamic percolation' algorithm. It's been on our todo list forever but it seems like a large undertaking and we are a bit selfish because we don't want to invest the time on it.