Chaning the Drainage to Imbibition process in the Predicting relative permeability

[JMCarmona10]

I'm encountering an issue with my code related to relative permeability prediction. The example code I'm working with simulates a drainage process, where a non-wetting phase (air) displaces a wetting phase (water). However, I want to modify this to simulate an imbibition process instead, where the wetting phase displaces the non-wetting phase. Additionally, I want to calculate the wetting phase saturation rather than the non-wetting phase saturation.

I've attempted to modify the code to simulate imbibition by adjusting the contact angle, but the resulting graph remains unchanged, which suggests there might be an error in my implementation.

Could you please help me identify what I'm doing wrong? I'll share the code for reference.

Relative Permeability Calculation using OpenPNM

This code simulates fluid invasion and calculates relative permeability using OpenPNM (Open Pore Network Modeling).

# Initialize the pore network
pn = op.network.Cubic(shape=[15, 15, 15], spacing=1e-6)
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()

# Set up air phase
air = op.phase.Air(network=pn, name='air')
air.add_model_collection(op.models.collections.phase.air)
air.add_model_collection(op.models.collections.physics.basic)
air.regenerate_models()

# Set up water phase with specific properties
water = op.phase.Water(network=pn, name='water')
water['pore.surface_tension'] = 15.72
water['pore.contact_angle'] = 30.0
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.basic)
water.regenerate_models()

# Apply invasion percolation
ip = op.algorithms.InvasionPercolation(network=pn, phase=water)
Finlets_init = pn.pores('left')
Finlets = ([Finlets_init[x] for x in range(0, len(Finlets_init), 2)])
ip.set_inlet_BC(pores=Finlets)
ip.run()

def sat_occ_update(network, nwp, wp, ip, i):
    """
    Calculates the saturation of wetting phase (S_wp) using the invasion
    sequence from invasion percolation.
    
    Parameters
    ----------
    network: network
    nwp : phase
        non-wetting phase
    wp : phase
        wetting phase
    ip : IP
        invasion percolation (ran before calling this function)
    i: int
        The invasion_sequence limit for masking pores/throats 
    """
    # Create masks for invaded pores/throats
    pore_mask = ip['pore.invasion_sequence'] < i
    throat_mask = ip['throat.invasion_sequence'] < i
    
    # Calculate volumes for wetting phase (non-invaded volumes)
    wp_sat_p = np.sum(network['pore.volume'][~pore_mask])
    wp_sat_t = np.sum(network['throat.volume'][~throat_mask])
    wp_sat = wp_sat_p + wp_sat_t
    bulk = network['pore.volume'].sum() + network['throat.volume'].sum()
    
    sat = wp_sat/bulk  # This calculates S_wp
    
    # Update phase occupancy
    nwp['pore.occupancy'] = 1-pore_mask
    nwp['throat.occupancy'] = 1-throat_mask
    wp['throat.occupancy'] = throat_mask
    wp['pore.occupancy'] = pore_mask
    
    return sat

def Rate_calc(network, phase, inlet, outlet, conductance):
    """
    Simulates stokes flow given the hydraulic conductance keyword 
    and returns the rate of inlet fluid flow.
    
    Parameters
    ----------
    network : network
    phase : phase
    inlet : array
        inlet pores
    outlet : array
        outlet pores
    conductance : str
        name of conductance model to use
    """
    phase.regenerate_models()
    St_p = op.algorithms.StokesFlow(network=network, phase=phase)
    St_p.settings._update({'conductance': conductance})
    St_p.set_value_BC(pores=inlet, values=1)
    St_p.set_value_BC(pores=outlet, values=0)
    St_p.run()
    val = np.abs(St_p.rate(pores=inlet, mode='group'))
    return val

# Set up flow parameters
flow_in = pn.pores('left')
flow_out = pn.pores('right')

# Define multiphase conductance
model_mp_cond = op.models.physics.multiphase.conduit_conductance
air.add_model(model=model_mp_cond, propname='throat.conduit_hydraulic_conductance',
              throat_conductance='throat.hydraulic_conductance', mode='medium', regen_mode='deferred')
water.add_model(model=model_mp_cond, propname='throat.conduit_hydraulic_conductance',
              throat_conductance='throat.hydraulic_conductance', mode='medium', regen_mode='deferred')

# Calculate relative permeabilities
Swp_num = 50
flow_in = pn.pores('left')
flow_out = pn.pores('right')
max_seq = np.max([np.max(ip['pore.invasion_sequence']),
                  np.max(ip['throat.invasion_sequence'])])
start = max_seq//Swp_num
stop = max_seq
step = max_seq//Swp_num

Swparr = []
relperm_nwp = []
relperm_wp = []

for i in range(start, stop, step):
    air.regenerate_models()
    water.regenerate_models()
    
    sat = sat_occ_update(network=pn, nwp=water, wp=air, ip=ip, i=i)
    Swparr.append(sat)
    
    Rate_abs_nwp = Rate_calc(pn, air, flow_in, flow_out, 
                            conductance='throat.hydraulic_conductance')
    Rate_abs_wp = Rate_calc(pn, water, flow_in, flow_out,
                           conductance='throat.hydraulic_conductance')
    Rate_enwp = Rate_calc(pn, air, flow_in, flow_out,
                         conductance='throat.conduit_hydraulic_conductance')
    Rate_ewp = Rate_calc(pn, water, flow_in, flow_out, 
                        conductance='throat.conduit_hydraulic_conductance')
    
    relperm_nwp.append(Rate_enwp/Rate_abs_nwp)
    relperm_wp.append(Rate_ewp/Rate_abs_wp)

# Plot the results
plt.figure(figsize=[6,5])
plt.plot(Swparr, relperm_nwp, '*-', label='Kr_water')
plt.plot(Swparr, relperm_wp, 'o-', label='Kr_air')
plt.xlabel('Swp')
plt.ylabel('Kr')
plt.title('Relative Permeability in x direction')
plt.legend()
plt.show()

Here are the results where im changing the contact angle and im getting the same results which is not correct.