Meniscus Model Comparison#

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:

https://media.giphy.com/media/AIbz7mpqxgc5a/giphy.gif

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.

https://i.imgur.com/b2der2w.png

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.

[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)
[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.

[3]:
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import openpnm.models.physics as pm
net = op.network.Cubic(shape=[2, 1, 1], spacing=5e-5)
geo = op.geometry.SpheresAndCylinders(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
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [3], in <cell line: 8>()
      4 net = op.network.Cubic(shape=[2, 1, 1], spacing=5e-5)
      5 geo = op.geometry.SpheresAndCylinders(network=net,
      6                                pores=net.pores(),
      7                                throats=net.throats())
----> 8 phase = op.phases.Water(network=net)
      9 phase['pore.contact_angle'] = theta
     10 phys = op.physics.Standard(network=net,
     11                            phase=phase,
     12                            geometry=geo)

AttributeError: module 'openpnm' has no attribute 'phases'

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.

[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)#

[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))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [5], in <cell line: 3>()
      1 circular_model = pm.meniscus.purcell
----> 3 phys.add_model(propname='throat.max',
      4                model=circular_model,
      5                mode='max',
      6                r_toroid=fiberRad)
      7 phys.add_model(propname='throat.touch',
      8                model=circular_model,
      9                mode='touch',
     10                r_toroid=fiberRad)
     11 phys.add_model(propname='throat.meniscus',
     12                model=circular_model,
     13                mode='men',
     14                r_toroid=fiberRad,
     15                target_Pc=1000)

NameError: name 'phys' is not defined

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.

[6]:
ax = plot_meniscus(target_Pc=touch_Pc, meniscus_model=circular_model)
ax.plot([max_bulge, max_bulge], [-throatRad, throatRad], 'r--');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 ax = plot_meniscus(target_Pc=touch_Pc, meniscus_model=circular_model)
      2 ax.plot([max_bulge, max_bulge], [-throatRad, throatRad], 'r--')

NameError: name 'touch_Pc' is not defined
[7]:
ax = plot_meniscus(target_Pc=max_Pc_circle, meniscus_model=circular_model)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 ax = plot_meniscus(target_Pc=max_Pc_circle, meniscus_model=circular_model)

NameError: name 'max_Pc_circle' is not defined

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

[8]:
interact(plot_meniscus,
         target_Pc=(-2000, max_Pc_circle, 1),
         meniscus_model=fixed(circular_model),
         ax=fixed(None));
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 1>()
      1 interact(plot_meniscus,
----> 2          target_Pc=(-2000, max_Pc_circle, 1),
      3          meniscus_model=fixed(circular_model),
      4          ax=fixed(None))

NameError: name 'max_Pc_circle' is not defined

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:

[9]:
sinusoidal_model = pm.meniscus.sinusoidal
[10]:
display(sinusoidal_model)
<function openpnm.models.physics.meniscus._funcs.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')>
[11]:
phys.add_model(propname='throat.meniscus',
               model=sinusoidal_model,
               mode='men',
               r_toroid=fiberRad,
               target_Pc=1000)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 phys.add_model(propname='throat.meniscus',
      2                model=sinusoidal_model,
      3                mode='men',
      4                r_toroid=fiberRad,
      5                target_Pc=1000)

NameError: name 'phys' is not defined

The equation for the solid sinusoidal profile is:

[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
[12]:
$\displaystyle - b \cos{\left(\frac{\pi x}{2 a} \right)} + b + rt$
[13]:
# Derivative of profile
rprime = r.diff(x)
rprime
[13]:
$\displaystyle \frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a}$
[14]:
# Filling angle
alpha = sym_atan(rprime)
alpha
[14]:
$\displaystyle \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)}$
[15]:
# Angle between y axis, meniscus center and meniscus contact point
eta = sym_pi - (theta + alpha)
eta
[15]:
$\displaystyle - \theta - \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} + \pi$
[16]:
# Angle between x axis, meniscus center and meniscus contact point
gamma = sym_pi/2 - eta
gamma
[16]:
$\displaystyle \theta + \operatorname{atan}{\left(\frac{\pi b \sin{\left(\frac{\pi x}{2 a} \right)}}{2 a} \right)} - \frac{\pi}{2}$
[17]:
# Radius of curvature of meniscus
rm = r/sym_cos(eta)
rm
[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)}}$
[18]:
# Distance along x-axis from center of curvature to meniscus contact point
d = rm*sym_sin(eta)
d
[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)}}$
[19]:
# Capillary Pressure
p = 2*sigma/rm
p
[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}$
[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)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [20], in <cell line: 1>()
----> 1 phys.add_model(propname='throat.max',
      2                model=sinusoidal_model,
      3                mode='max',
      4                r_toroid=fiberRad)
      5 phys.add_model(propname='throat.touch',
      6                model=sinusoidal_model,
      7                mode='touch',
      8                r_toroid=fiberRad)
      9 max_Pc_sin = phys['throat.max'][0]

NameError: name 'phys' is not defined
[21]:
plot_meniscus(target_Pc=max_Pc_sin, meniscus_model=sinusoidal_model);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [21], in <cell line: 1>()
----> 1 plot_meniscus(target_Pc=max_Pc_sin, meniscus_model=sinusoidal_model)

NameError: name 'max_Pc_sin' is not defined
[22]:
interact(plot_meniscus,
         target_Pc=(-2000, max_Pc_sin, 1),
         meniscus_model=fixed(sinusoidal_model),
         ax=fixed(None));
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [22], in <cell line: 1>()
      1 interact(plot_meniscus,
----> 2          target_Pc=(-2000, max_Pc_sin, 1),
      3          meniscus_model=fixed(sinusoidal_model),
      4          ax=fixed(None))

NameError: name 'max_Pc_sin' is not defined

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:

[23]:
y = sym_sqrt(1 - (x/a)**2)*b
y
[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

[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)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [24], in <cell line: 2>()
      1 # Scale ellipse in x direction
----> 2 phys['throat.scale_a'] = fiberRad
      3 # Scale ellipse in y direction
      4 phys['throat.scale_b'] = fiberRad

NameError: name 'phys' is not defined
[25]:
plot_meniscus(target_Pc=max_Pc_ellipse, meniscus_model=general_model);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [25], in <cell line: 1>()
----> 1 plot_meniscus(target_Pc=max_Pc_ellipse, meniscus_model=general_model)

NameError: name 'max_Pc_ellipse' is not defined
[26]:
max_Pc_ellipse
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [26], in <cell line: 1>()
----> 1 max_Pc_ellipse

NameError: name 'max_Pc_ellipse' is not defined
[27]:
interact(plot_meniscus,
         target_Pc=(-2000, max_Pc_ellipse, 1),
         meniscus_model=fixed(general_model),
         ax=fixed(None));
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [27], in <cell line: 1>()
      1 interact(plot_meniscus,
----> 2          target_Pc=(-2000, max_Pc_ellipse, 1),
      3          meniscus_model=fixed(general_model),
      4          ax=fixed(None))

NameError: name 'max_Pc_ellipse' is not defined

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.

[28]:
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])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [28], in <cell line: 2>()
      1 bs = np.linspace(0.2, 1.0, 4)*throatRad
----> 2 phys['throat.scale_a'] = throatRad
      3 elliptical_pressures = []
      4 sinusoidal_pressures = []

NameError: name 'phys' is not defined
[29]:
plt.figure()
plt.plot(bs/throatRad, elliptical_pressures, 'g-');
plt.plot(bs/throatRad, sinusoidal_pressures, 'r-');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [29], in <cell line: 2>()
      1 plt.figure()
----> 2 plt.plot(bs/throatRad, elliptical_pressures, 'g-');
      3 plt.plot(bs/throatRad, sinusoidal_pressures, 'r-')

NameError: name 'elliptical_pressures' is not defined
<Figure size 360x360 with 0 Axes>

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.