Part D: Comparison of toroidal meniscus models with different profile shapes

Introduction

So far all the capillary entry pressures for the percoaltion examples were calculated using the Standard physics model which is the Washburn model for straight walled capillary tubes. This has been shown to be a bad model for fibrous media where the walls of throats are converging and diverging. In the study Capillary Hysteresis in Neutrally Wettable Fibrous Media: A Pore Network Study of a Fuel Cell Electrode percolation in fibrous media was simulated using a meniscus model that assumed the contrictions between fibers are similar to a toroid:

This model was first proposed by Purcell and treats the inner solid profile as a circle. As the fluid invades through the center of the torus the meniscus is pinned to the surface and the "effective" contact angle becomes influenced by the converging diverging geometry and is a function of the filling angle $\alpha$. The shape of the meniscus as the invading phase moves upwards through the torus with key model parameters is shown below.

Different intrinsic contact angles through invading phase are shown above: (a) 60$^\circ$, (b) 90$^\circ$ and (c) 120$^\circ$. All scenarios clearly show an inflection of the meniscus curvature signifying a switch in the sign of the capillary pressure from negative to positive. This inflection is predicted to occur for all contact angles by the model with varying filling angle. The capillary pressure can be shown to be:

$P_C = -2\sigma cos(\theta-\alpha))/(r+R(1-cos(\alpha))$

A consequence of the circular solid profile is that all fluid behaves as non-wetting fluid because $\alpha$ can range from -90$^\circ$ to 90$^\circ$ degrees and so even if $\theta$ is 0 then the meniscus is still pinned at zero capillary pressure at the very furthest part of the throat where the $\alpha$ is 90$^\circ$

Considering other shapes of solid profile this situation can be avoided. It will be shown by reformulating the Purcell model in a more general way that allows for a flexible defintion of the solid profile that filling angle can be limited to values below 90 and allow for spontaneous imbibition (percolation threshold below zero) of highly wetting fluids.

Set up

We will set up a trivially small network with one throat to demonstrate the use of the meniscus model. Here we do the imports and define a few functions for plotting.



In [1]:

    
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import sympy as syp
from sympy import lambdify, symbols
from sympy import atan as sym_atan
from sympy import cos as sym_cos
from sympy import sin as sym_sin
from sympy import sqrt as sym_sqrt
from sympy import pi as sym_pi
from ipywidgets import interact, fixed
from IPython.display import display
import warnings
np.random.seed(10)
warnings.simplefilter(action='ignore')
matplotlib.rcParams['figure.figsize'] = (5, 5)



In [2]:

    
theta = 60
fiberRad = 5e-6
throatRad = 1e-5
max_bulge = 1e-5

Now we define our two pore network and add the meniscus model in several modes: 'max' returns the maximum pressure experienced by the meniscus as it transitions through the throat, i.e. the burst entry pressure. 'touch' is the pressure at which the meniscus has protruded past the throat center a distance defined by the 'touch_length' dictionary key. In network simulations this could be set to the pore_diameter. Finally the 'men' mode accepts a target_Pc parameter and returns all the mensicus information required for assessing cooperative filling or plotting.



In [3]:

    
import openpnm as op
import openpnm.models.physics as pm
net = op.network.Cubic(shape=[2, 1, 1], spacing=5e-5)
geo = op.geometry.StickAndBall(network=net,
                               pores=net.pores(),
                               throats=net.throats())
phase = op.phases.Water(network=net)
phase['pore.contact_angle'] = theta
phys = op.physics.Standard(network=net,
                           phase=phase,
                           geometry=geo)
geo['throat.diameter'] = throatRad*2
geo['throat.touch_length'] = max_bulge

We define a plotting function that uses the meniscus data: $\alpha$ is filling angle as defined above, $radius$ is the radius of curvature of the mensicus, $center$ is the position of the centre of curvature relative to the throat center along the axis of the throat, $\gamma$ is the angle between the throat axis and the line joining the meniscus center and meniscus contact point.



In [4]:

    
def plot_meniscus(target_Pc, meniscus_model=None, ax=None):
    throatRad = geo['throat.diameter'][0]/2
    theta = np.deg2rad(phys['pore.contact_angle'][0])
    throat_a = phys['throat.scale_a']
    throat_b = phys['throat.scale_b']
    x_points = np.arange(-0.99, 0.99, 0.01)*throat_a
    if ax is None:
        fig, ax = plt.subplots()
    if meniscus_model.__name__ == 'purcell':
        # Parameters for plotting fibers
        x, R, rt, s, t = syp.symbols('x, R, rt, s, t')
        y = R*syp.sqrt(1- (x/R)**2)
        r = rt + (R-y)
        rx = syp.lambdify((x, R, rt), r, 'numpy')
        ax.plot(x_points, rx(x_points, fiberRad, throatRad), 'k-');
        ax.plot(x_points, -rx(x_points, fiberRad, throatRad), 'k-');
        phys.add_model(propname='throat.meniscus',
                       model=meniscus_model,
                       mode='men',
                       r_toroid=fiberRad,
                       target_Pc=target_Pc)
    elif meniscus_model.__name__ == 'sinusoidal':
        x, a, b, rt, sigma, theta = syp.symbols('x, a, b, rt, sigma, theta')
        y = (sym_cos(sym_pi*x/(2*a)))*b
        r = rt + (b-y)
        rx = lambdify((x, a, b, rt), r, 'numpy')
        ax.plot(x_points, rx(x_points, throat_a, throat_b, throatRad), 'k-');
        ax.plot(x_points, -rx(x_points, throat_a, throat_b, throatRad), 'k-');
        phys.add_model(propname='throat.meniscus',
                       model=meniscus_model,
                       mode='men',
                       r_toroid=fiberRad,
                       target_Pc=target_Pc)
    else:
        # General Ellipse
        x, a, b, rt, sigma, theta = syp.symbols('x, a, b, rt, sigma, theta')
        profile_equation = phys.models['throat.entry_pressure']['profile_equation']
        if profile_equation == 'elliptical':
            y = sym_sqrt(1 - (x/a)**2)*b
        elif profile_equation == 'sinusoidal':
            y = (sym_cos(sym_pi*x/(2*a)))*b
        r = rt + (b-y)
        rx = lambdify((x, a, b, rt), r, 'numpy')
        ax.plot(x_points, rx(x_points, throat_a, throat_b, throatRad), 'k-');
        ax.plot(x_points, -rx(x_points, throat_a, throat_b, throatRad), 'k-');
        phys.add_model(propname='throat.meniscus',
                       model=meniscus_model,
                       profile_equation=profile_equation,
                       mode='men',
                       target_Pc=target_Pc)
    men_data = {}
    men_data['alpha'] = phys['throat.meniscus.alpha']
    men_data['gamma'] = phys['throat.meniscus.gamma']
    men_data['radius'] = phys['throat.meniscus.radius']
    men_data['center'] = phys['throat.meniscus.center']
    arc_cen = men_data['center']
    arc_rad = men_data['radius']
    arc_angle = men_data['gamma']
    angles = np.linspace(-arc_angle, arc_angle, 100)
    arcx = arc_cen + arc_rad*np.cos(angles)
    arcy = arc_rad*np.sin(angles)
    ax.plot(arcx, arcy, 'b-')
    ax.scatter(phys['throat.meniscus.pos'], phys['throat.meniscus.rx']);
    ax.axis('equal')
    ax.ticklabel_format(style='sci', axis='both', scilimits=(-6,-6))
    return ax

Circular (Purcell)



In [5]:

    
circular_model = pm.meniscus.purcell

phys.add_model(propname='throat.max',
               model=circular_model,
               mode='max',
               r_toroid=fiberRad)
phys.add_model(propname='throat.touch',
               model=circular_model,
               mode='touch',
               r_toroid=fiberRad)
phys.add_model(propname='throat.meniscus',
               model=circular_model,
               mode='men',
               r_toroid=fiberRad,
               target_Pc=1000)
touch_Pc = phys['throat.touch'][0]
print('Pressure at maximum bulge', np.around(touch_Pc, 0))
max_Pc_circle = phys['throat.max'][0]
print('Circular profile critical entry pressure', np.around(max_Pc_circle, 0))









    



Pressure at maximum bulge 7213.0
Circular profile critical entry pressure 8165.0

We can see that the touch_Pc calculated earlier, corresponds with the tip of the meniscus exceeding the max_bulge parameter. Try changing this and re-running to see what happens.



In [6]:

    
#NBVAL_IGNORE_OUTPUT
ax = plot_meniscus(target_Pc=touch_Pc, meniscus_model=circular_model)
ax.plot([max_bulge, max_bulge], [-throatRad, throatRad], 'r--');



In [7]:

    
#NBVAL_IGNORE_OUTPUT
ax = plot_meniscus(target_Pc=max_Pc_circle, meniscus_model=circular_model)

We can interact with the mensicus model by changing the target_Pc parameter.



In [8]:

    
#NBVAL_IGNORE_OUTPUT
interact(plot_meniscus,
         target_Pc=(-2000, max_Pc_circle, 1),
         meniscus_model=fixed(circular_model),
         ax=fixed(None));

Here we can see that the critical entry pressure for the circular profile is positive, even though the intrinsic contact angle is highly non-wetting

Sinusoidal

Now we can start to compare the different meniscus models:



In [9]:

    
sinusoidal_model = pm.meniscus.sinusoidal



In [10]:

    
display(sinusoidal_model)









    





<function openpnm.models.physics.meniscus.sinusoidal(target, mode='max', target_Pc=None, num_points=1000.0, r_toroid=5e-06, throat_diameter='throat.diameter', pore_diameter='pore.diameter', touch_length='throat.touch_length', surface_tension='pore.surface_tension', contact_angle='pore.contact_angle')>



In [11]:

    
phys.add_model(propname='throat.meniscus',
               model=sinusoidal_model,
               mode='men',
               r_toroid=fiberRad,
               target_Pc=1000)

The equation for the solid sinusoidal profile is:



In [12]:

    
x, a, b, rt, sigma, theta = syp.symbols('x, a, b, rt, sigma, theta')
y = (sym_cos(sym_pi*x/(2*a)))*b
r = rt + b-y
r









    Out[12]:





$\displaystyle - b \cos{\left(\frac{\pi x}{2 a} \right)} + b + rt$



In [13]:

    
# Derivative of profile
rprime = r.diff(x)
rprime









    Out[13]:





$\displaystyle \frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a}$



In [14]:

    
# Filling angle
alpha = sym_atan(rprime)
alpha









    Out[14]:





$\displaystyle \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)}$



In [15]:

    
# angle between y axis, meniscus center and meniscus contact point
eta = sym_pi - (theta + alpha)
eta









    Out[15]:





$\displaystyle - \theta - \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} + \pi$



In [16]:

    
# angle between x axis, meniscus center and meniscus contact point
gamma = sym_pi/2 - eta
gamma









    Out[16]:





$\displaystyle \theta + \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} - \frac{\pi}{2}$



In [17]:

    
# Radius of curvature of meniscus
rm = r/sym_cos(eta)
rm









    Out[17]:





$\displaystyle - \frac{- b \cos{\left(\frac{\pi x}{2 a} \right)} + b + rt}{\cos{\left(\theta + \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} \right)}}$



In [18]:

    
# distance along x-axis from center of curvature to meniscus contact point
d = rm*sym_sin(eta)
d









    Out[18]:





$\displaystyle - \frac{\left(- b \cos{\left(\frac{\pi x}{2 a} \right)} + b + rt\right) \sin{\left(\theta + \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} \right)}}{\cos{\left(\theta + \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} \right)}}$



In [19]:

    
# Capillary Pressure
p = 2*sigma/rm
p









    Out[19]:





$\displaystyle - \frac{2 \sigma \cos{\left(\theta + \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} \right)}}{- b \cos{\left(\frac{\pi x}{2 a} \right)} + b + rt}$



In [20]:

    
phys.add_model(propname='throat.max',
               model=sinusoidal_model,
               mode='max',
               r_toroid=fiberRad)
phys.add_model(propname='throat.touch',
               model=sinusoidal_model,
               mode='touch',
               r_toroid=fiberRad)
max_Pc_sin = phys['throat.max'][0]
print(max_Pc_sin)









    



4729.770413396985



In [21]:

    
#NBVAL_IGNORE_OUTPUT
plot_meniscus(target_Pc=max_Pc_sin, meniscus_model=sinusoidal_model);



In [22]:

    
#NBVAL_IGNORE_OUTPUT
interact(plot_meniscus,
         target_Pc=(-2000, max_Pc_sin, 1),
         meniscus_model=fixed(sinusoidal_model),
         ax=fixed(None));

Now the crtical entry pressure is negative signifying that spontaneous imbibition will occur

General Elliptical

Similarly we can define an elliptical profile and use the same method to determine the capillary pressure:



In [23]:

    
y = sym_sqrt(1 - (x/a)**2)*b
y









    Out[23]:





$\displaystyle b \sqrt{1 - \frac{x^{2}}{a^{2}}}$

In-fact this is the model that OpenPNM uses for Purcell as well with a = b = fiber radius



In [24]:

    
# Scale ellipse in x direction
phys['throat.scale_a'] = fiberRad
# Scale ellipse in y direction
phys['throat.scale_b'] = fiberRad
general_model = pm.meniscus.general_toroidal
phys.add_model(propname='throat.entry_pressure',
               model=general_model,
               profile_equation='elliptical',
               mode='max')
max_Pc_ellipse = phys['throat.entry_pressure'][0]
print(max_Pc_ellipse)









    



8165.324889242946



In [25]:

    
#NBVAL_IGNORE_OUTPUT
plot_meniscus(target_Pc=max_Pc_ellipse, meniscus_model=general_model);



In [26]:

    
max_Pc_ellipse









    Out[26]:





8165.324889242946



In [27]:

    
#NBVAL_IGNORE_OUTPUT
interact(plot_meniscus,
         target_Pc=(-2000, max_Pc_ellipse, 1),
         meniscus_model=fixed(general_model),
         ax=fixed(None));

The two scale factors can now be used to determine a wide range of capillary behaviours with one general model. Below we run the model for a range of scaling factors showing the effect on the sign and magnitude of the entry pressure.



In [28]:

    
#NBVAL_IGNORE_OUTPUT
bs = np.linspace(0.2, 1.0, 4)*throatRad
phys['throat.scale_a'] = throatRad
elliptical_pressures = []
sinusoidal_pressures = []
fig, (ax1, ax2) = plt.subplots(2, len(bs), figsize=(10, 10))
for i in range(len(bs)):
    phys['throat.scale_b'] = bs[i]
    phys.add_model(propname='throat.entry_pressure',
                   model=general_model,
                   profile_equation='elliptical',
                   mode='max',
                   num_points=1000)
    Pc = phys['throat.entry_pressure']
    elliptical_pressures.append(Pc)
    plot_meniscus(target_Pc=Pc, meniscus_model=general_model, ax=ax1[i])
for i in range(len(bs)):
    phys['throat.scale_b'] = bs[i]
    phys.add_model(propname='throat.entry_pressure',
                   model=general_model,
                   profile_equation='sinusoidal',
                   mode='max',
                   num_points=1000)
    Pc = phys['throat.entry_pressure']
    sinusoidal_pressures.append(Pc)
    plot_meniscus(target_Pc=Pc, meniscus_model=general_model, ax=ax2[i])



In [29]:

    
#NBVAL_IGNORE_OUTPUT
plt.figure()
plt.plot(bs/throatRad, elliptical_pressures, 'g-');
plt.plot(bs/throatRad, sinusoidal_pressures, 'r-');

Here we can see that the two different shaped profiles lead to quite different capiallary behaviour. The elliptical profile always resuls in positive pressure and the meniscus is basically pinned to the end of the throat where highest pressure occurs as alpha always reaches 90. Whereas the sinusiodal model allows for spontaneous imbibition where a breakthrough may occur at negative capillary pressure for wetting fluids if the wall angle is shallow.