[OpenPNM] Part 2 ——Geometric Properties

发布于 2023-03-30  518 次阅读


Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

指定几何属性

几何属性是指孔和流道的物理尺寸,例如孔径和流道长度。这些值对于网络建模至关重要,因为它们控制着网络的所有传输和渗透行为。与可以归一化的相属性不同,几何属性比连通性等拓扑属性更能定义网络。

几何属性可以通过各种方法计算,为孔隙分配随机值,从断层扫描图像中提取值等等。OpenPNM提供了一个函数库,或“孔隙尺度模型”,可以自动计算这些几何属性。本教程将涵盖以下主题:

  • 手动计算孔和流道的几何属性(不推荐,太麻烦)

  • 使用库中的孔隙尺度模型(推荐)

  • 依赖项处理程序概述

先创建一个简单的网络模型(只包含孔和流道的位置、数量信息,不包含几何参数)

import numpy as np
import matplotlib.pyplot as plt
import openpnm as op
op.visualization.set_mpl_style()
np.random.seed(0)
pn = op.network.Cubic([20, 20, 20], spacing=5e-5)
print(pn)

想要为网络模型添加几何参数可以通过下面几种方式:

手动计算几何属性

注意,手动计算非常麻烦,但是能很好的说明了计算过程,不推荐在编码过程中采用手动计算

让我们从添加孔和流道直径分布开始。有几种不同的方法可以做到这一点,我们将探索每种方法。

从分布中添加孔流道直径

scipy.stats模块定义了许多统计学分布,首先,让我们使用正态分布来生成1到50 um范围内的孔径

np.random.seed(0)
import scipy.stats as spst
psd = spst.norm.rvs(loc=25, scale=6, size=pn.Np)
plt.hist(psd, edgecolor='k')
plt.xlim([0, 50]);

Definition : rvs(*args, **kwds)
Random variates of given type.

Parameters

arg1, arg2, arg3,...array_like
The shape parameter(s) for the distribution (see docstring of the instance object for more information).

loc:array_like, optional
Location parameter (default=0).

scale:array_like, optional
Scale parameter (default=1).

size:int or tuple of ints, optional
Defining number of random variates (default is 1).

random_state:{None, int, numpy.random.Generator,
numpy.random.RandomState}, optional

If seed is None (or np.random), the numpy.random.RandomState singleton is used. If seed is an int, a new RandomState instance is used, seeded with seed. If seed is already a Generator or RandomState instance then that instance is used.

Returns

rvsndarray or scalar
Random variates of given size.

plt.hist()相关说明可以参考:https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.hist.html
和https://blog.csdn.net/weixin_45520028/article/details/113924866

#上面的分布看起来不错,但是要确保我们的分布不是太宽,以至于它的值大于 50 um:

print(psd.min())
print(psd.max())
'''
2.5593961722893255
47.80996128980269
'''
#现在,我们将这些值转换单位并分配给网络:

pn['pore.diameter'] = psd*1e-6  # um to m

接下来定义流道直径。我们可以用同样的方式做到这一点,但使用不同统计学分布。请注意,不建议使用此方法,因为正如我们将看到的,它会导致流道直径大于它们连接两个的两个孔。我们将在后面解决此问题:

np.random.seed(0)
tsd = spst.weibull_min.rvs(c=1.5, loc=.5, scale=7.5, size=pn.Nt)
plt.hist(tsd, edgecolor='k')
plt.xlim([0, 50]);

spst.weibull_min.rvs()函数的相关信息可以参考:https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.weibull_min.html
和https://blog.csdn.net/weixin_41102672/article/details/106565711

#同样,让我们检查高值和低值:

print(tsd.min())
print(tsd.max())
'''
0.5130345457395142
36.96861960231873
'''
#所以我们可以看到我们的流道直径小到500纳米,大到37微米。这些也可以分配给网络:

pn['throat.diameter'] = tsd*1e-6  # um to m

这种方法的问题在于,孔和流道直径的大小只是分配给随机位置,而不是将小流道放在小孔之间,反之亦然。当流道直径大于它们连接的孔时,可能会导致结果出现问题,因为我们通常认为流道直径比两个孔之都小。可以通过下面的代码数一数有多少是有问题的

hits = np.any(pn['pore.diameter'][pn.conns].T < pn['throat.diameter'], axis=0)
hits.sum()
'''
564
'''

数组后加 .T 意思是矩阵转置

通过conns索引孔隙属性:使用数组索引孔隙属性将返回一个Nt-by-2数组,其中包含每个列中流道末端的孔隙属性。例如,看看流道两端孔的直径。connspn['pore.diameter'][pn.conns]

np.any()函数的相关信息可以参考:https://numpy.org/doc/stable/reference/generated/numpy.any.html

有两种更好的方法来分配孔和流道的直径,以确保流道直径小于它们连接的孔径。第一种是通过随机分配孔径进行上述操作,然后将流道分配为与其连接的最小孔一样小(或更小)。这可以使用一些函数来完成,如下所示:numpy

tsd = np.amin(pn['pore.diameter'][pn.conns], axis=1)
plt.hist(tsd, edgecolor='k', density=True, alpha=0.5)#alpha设置透明度,0为完全透明
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.xlim([0, 50e-6]);

从上图中我们可以看到,流道直径大小与孔径相似,但是偏小一点这是意料之中的,因为它们取自孔径值,而且流道直径要比孔径小。

还需要通过某种因素迫使喉咙小于相邻的两个毛孔:

pn['throat.diameter'] = tsd*0.5
plt.hist(pn['throat.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.xlim([0, 50e-6]);

上述方法的一个缺点是无法再控制流道直径大小的分布。这可以通过应用于随机种子来解决。将随机种子放入累积密度函数中以计算相应的直径大小。如下所示:

随机种子的相关解释可以参考:https://blog.csdn.net/ding_programmer/article/details/95097924
https://en.wikipedia.org/wiki/Random_seed
https://www.geeksforgeeks.org/random-seed-in-python/

lo, hi = 0.01, 0.99
pn['pore.seed'] = np.random.rand(pn.Np)*(hi - lo) + lo
#设置种子的和值可以防止位于分布尾部较远的值出现,从而导致过大的孔径或喉道直径

#使用点概率函数返回给定概率下的累积分布值,因此上面计算的种子值可以使用以下代码与大小值相关联
import scipy.stats as spst
psd = spst.norm.ppf(loc=25, scale=6, q=pn['pore.seed'])
pn['pore.diameter'] = psd*1e-6
plt.hist(pn['pore.diameter'], edgecolor='k')
plt.xlim([0, 50e-6]);

Definition : ppf(q, *args, **kwds)
Percent point function (inverse of cdf) at q of the given RV.

Parameters

q:array_like
lower tail probability

arg1, arg2, arg3,...array_like
The shape parameter(s) for the distribution (see docstring of the instance object for more information)

loc:array_like, optional
location parameter (default=0)

scale:array_like, optional
scale parameter (default=1)

Returns

x:array_like
quantile corresponding to the lower tail probability q.

流道种子值为两个相邻孔隙的最小孔隙种子值:

pn['throat.seed'] = np.amin(pn['pore.seed'][pn.conns], axis=1)

同样,可以根据流道的种子算出流道的直径大小:

tsd = spst.norm.ppf(loc=25, scale=6, q=pn['throat.seed'])
pn['throat.diameter'] = tsd*1e-6
plt.hist(pn['throat.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.xlim([0, 50e-6]);

指定孔径和流道直径后,下一步是定义任何相关属性,例如孔体积或流道长度。这些额外的计算往往要求我们假设孔和流道的形状。通常是球棍模型。但是我们还可以使用立方体孔和长方体喉咙。出于探讨的原因,立方体形状的选择也恰好使计算变得更加容易,因此在本教程中,我们将假设立方体形状。让我们计算流道长度和孔表面积,因为它们中的每一个都说明了有关使用矢量化计算的重要概念。

计算流道长度

流道长是孔间间距,减去每个孔的半径。可以使用以下矢量化技巧在不需要 for 循环的情况下找到这一点

R1, R2 = (pn['pore.diameter'][pn.conns]/2).T
L_total = np.sqrt(np.sum(np.diff(pn.coords[pn.conns], axis=1).squeeze()**2, axis=1))
'''
np.diff():Calculate the n-th discrete difference along the given axis.
np.sum():Sum of array elements over a given axis.
numpy.squeeze() function is used when we want to remove single-dimensional entries from the shape of an array.
'''
Lt = L_total - R1 - R2
print(L_total)

上述单元格包含3个复杂但强大的步骤:

每个流道两端的孔半径被检索出来。因为每个孔都有很多流道,所以每个流道要一次孔半径。这在第3步中使用。R1R2

孔到孔的距离是通过检索每个孔的坐标来找到的,与计算半径的方法相同。然后计算每对孔隙之间的欧几里得距离,方法是找到它们的差值,对其求和,然后取平方根。

每个流道的长度是用孔之间的间距减去两个孔半径

这里的主要信息是流道属性可以用矢量化的方式计算,即使需要一些孔隙属性。

计算孔隙表面积

孔表面积是立方孔体的总面积,减去每个相邻流道的横截面积。这也可以在没有 for 循环的情况下找到,但需要了解一些更深层次的功能,

At = pn['throat.diameter']**2
SAp = (pn['pore.diameter']**2)*6
np.subtract.at(SAp, pn.conns[:, 0], At)
np.subtract.at(SAp, pn.conns[:, 1], At)
print(SAp)
'''
[[   0    1]
 [   1    2]
 [   2    3]
 ...
 [7597 7997]
 [7598 7998]
 [7599 7999]]
[7.27174552e-09 3.85799662e-09 4.84581355e-09 ... 2.94168620e-09
 2.75016941e-09 3.47031746e-09]
'''

上面的单元格包含2行重要性,每一行几乎都做同样的事情:

1、从数组的第一列中列出的每个孔中减去喉道横截面积。由于一个给定的孔隙与潜在的许多喉道相连,因此必须使用减法函数的特殊版本来完成减法函数,该函数从所指示的位置中减去每一个值,并且以累积的方式这样做。因此,如果孔隙i出现了不止一次,则删除的值就不止一个。该操作在原地执行,因此不会返回任何数组

2、从数组的第二列中列出的每个孔中减去喉道横截面积

使用OpenPNM自带的pore-scale模型

在本节中,我们将利用OpenPNM的“孔隙尺度”模型库来更轻松地完成相同的工作。首先生成一个新网络:

np.random.seed(0)
pn = op.network.Cubic([20, 20, 20], spacing=5e-5)

#为每个孔分配一个随机数,使用孔隙尺度模型:np.random.rand

f = op.models.geometry.pore_seed.random
pn.add_model(propname='pore.seed', 
             model=f,
             seed=None,  # Seed value provided to numpy's random number generator.
             num_range=[0.01, 0.99],)

#让我们让每个流道继承其每个相邻毛孔中的最小种子值,OpenPNM 为此提供了一个模型:
f = op.models.geometry.throat_seed.from_neighbor_pores
pn.add_model(propname='throat.seed', 
             model=f,
             prop='pore.seed',
             mode='min')
"""
add_model(propname, model, domain='all', regen_mode='normal', *, map: Mapping[_KT, _VT], iterable: Iterable[Tuple[_KT, _VT]], **kwargs: _VT)

        Add a pore-scale model to the object, along with the desired arguments

        Parameters
        ----------
        propname : str
            The name of the property being computed. E.g. if
            ``propname='pore.diameter'`` then the computed results will be stored
            in ``obj['pore.diameter']``.
        model : function handle
            The function that produces the values
        domain : str
            The label indicating the locations in which results generated by
            ``model`` should be stored. See `Notes` for more details.
        regen_mode : str
            How the model should be regenerated. Options are:

            ============ =====================================================
            regen_mode   description
            ============ =====================================================
            normal       (default) The model is run immediately upon being
                         added, and is also run each time ``regenerate_models``
                         is called.
            deferred     The model is NOT run when added, but is run each time
                         ``regenerate_models`` is called. This is useful for
                         models that depend on other data that may not exist
                         yet.
            constant     The model is run immediately upon being added, but is
                         is not run when ``regenerate_models`` is called,
                         effectively turning the property into a constant.
            ============ =====================================================

        kwargs : keyword arguments
            All additional keyword arguments are passed on to the model

        Notes
        -----
        The ``domain`` argument dictates where the results of ``model`` should
        be stored. For instance, given ``propname='pore.diameter'`` and
        ``domain='left'`` then when `model` is run, the results are stored in
        in the pores labelled left. Note that if ``model`` returns ``Np``
        values, then values not belonging to ``'pore.left'`` are discarded.
        The following steps outline the process:

        1. Find the pore indices:

        .. code-block:: python

          Ps = obj.pores('left')

        2. Run the model:

        .. code-block:: python

          vals = model(**kwargs)

        3. If the model returns a full Np-length array, then extract the
        correct values and apply them to the corresponding locations:

        .. code-block:: python

          if len(vals) == obj.Np:
              obj['pore.diameter'][Ps] = vals[Ps]

        4. If the model was designed to return only the subset of values then:

        .. code-block:: python

          if len(vals) == obj.num_pores('left'):
              obj['pore.diameter'][Ps] = vals

"""

现在我们可以从分布中计算孔和喉部大小,这些分布也以 OpenPNM 模型的形式提供:

f = op.models.geometry.pore_size.normal
pn.add_model(propname='pore.diameter',
             model=f,
             scale=1e-5, 
             loc=2.5e-5,
             seeds='pore.seed')
plt.hist(pn['pore.diameter'], edgecolor='k');

类似地,我们也获取流道大小的模型:

f = op.models.geometry.throat_size.normal
pn.add_model(propname='throat.diameter',
             model=f,
             scale=1e-5, 
             loc=2.5e-5)
plt.hist(pn['pore.diameter'], edgecolor='k', density=True, alpha=0.5)
plt.hist(pn['throat.diameter'], edgecolor='k', density=True, alpha=0.5);

现在我们有了孔和流道的直径的大小,我们可以计算它们的体积(和其他属性)

f = op.models.geometry.throat_length.spheres_and_cylinders
pn.add_model(propname='throat.length',
             model=f)
f1 = op.models.geometry.pore_volume.sphere 
pn.add_model(propname='pore.volume',
             model=f1)
f2 = op.models.geometry.throat_volume.cylinder
pn.add_model(propname='throat.volume',
             model=f2)

#球形孔与其连接的流道之间存在重叠区域。需要额外的步骤来消除:
f = op.models.geometry.throat_length.spheres_and_cylinders
pn.add_model(propname='throat.length',
             model=f)
f1 = op.models.geometry.pore_volume.sphere 
pn.add_model(propname='pore.volume',
             model=f1)
f2 = op.models.geometry.throat_volume.cylinder
pn.add_model(propname='throat.total_volume',
             model=f2)
f3 = op.models.geometry.throat_volume.lens
pn.add_model(propname='throat.lens_volume',
             model=f3)
f4 = op.models.misc.difference
pn.add_model(propname='throat.volume',
             model=f4,
             props=['throat.total_volume', 'throat.lens_volume'])
#检查一下
print("Volumes of full throats:", pn['throat.total_volume'][:3])
print("Volumes of lenses:", pn['throat.lens_volume'][:3])
print("Actual throat volumes:", pn['throat.volume'][:3])
'''
Volumes of full throats: [2.27203214e-14 2.58686910e-14 2.44028684e-14]
Volumes of lenses: [6.92245838e-15 8.40091924e-15 7.59556144e-15]
Actual throat volumes: [1.57978631e-14 1.74677718e-14 1.68073070e-14]
'''

使用预定义的模型集合

在 OpenPNM V3 中,我们引入了模型集合的概念 openpnm.models.collections,模型集合是形成完整且正确的几何形状的模型的预定义字典。例如,有一个名为 modelsspheres_and_cylinders 的集合,其中包含描述此几何图形所需的所有模型:

pn = op.network.Cubic(shape=[20, 20, 20], spacing=5e-5)
#The models collection is fetched as follows:

from pprint import pprint
mods = op.models.collections.geometry.spheres_and_cylinders
pprint(mods)
'''
{'pore.diameter': {'model': ,
                   'props': ['pore.max_size', 'pore.seed']},
 'pore.max_size': {'iters': 10,
                   'model': },
 'pore.seed': {'element': 'pore',
               'model': ,
               'num_range': [0.2, 0.7],
               'seed': None},
 'pore.volume': {'model': ,
                 'pore_diameter': 'pore.diameter'},
 'throat.cross_sectional_area': {'model': ,
                                 'throat_diameter': 'throat.diameter'},
 'throat.diameter': {'factor': 0.5,
                     'model': ,
                     'prop': 'throat.max_size'},
 'throat.diffusive_size_factors': {'model': ,
                                   'pore_diameter': 'pore.diameter',
                                   'throat_diameter': 'throat.diameter'},
 'throat.hydraulic_size_factors': {'model': ,
                                   'pore_diameter': 'pore.diameter',
                                   'throat_diameter': 'throat.diameter'},
 'throat.length': {'model': ,
                   'pore_diameter': 'pore.diameter',
                   'throat_diameter': 'throat.diameter'},
 'throat.lens_volume': {'model': ,
                        'pore_diameter': 'pore.diameter',
                        'throat_diameter': 'throat.diameter'},
 'throat.max_size': {'mode': 'min',
                     'model': ,
                     'prop': 'pore.diameter'},
 'throat.total_volume': {'model': ,
                         'throat_diameter': 'throat.diameter',
                         'throat_length': 'throat.length'},
 'throat.volume': {'model': ,
                   'props': ['throat.total_volume', 'throat.lens_volume']
'''

字典的打印输出不是很好,但可以看出键是模型名称,值是模型本身和模型的各种参数。可以使用 add_models_collection 将此字典添加到网络中,如下所示

pn.add_model_collection(mods)
pn.regenerate_models()
print(pn)

通过覆盖模型或其参数来自定义模型

OpenPNM 中包含的预先编写的模型和集合很有用,但通常需要调整某些模型以适应特定的应用程序。这可以通过更改现有模型的参数或将它们全部覆盖来完成。

让我们首先看看如何用我们选择的模型替换一个模型。该集合使用孔径的随机分布,以0和到最近孔的距离为边界。这可以通过打印具体的模型来检查,如下所示:
spheres_and_cylinders

print(pn.models['pore.diameter@all'])

用正态分布替换它:

f = op.models.geometry.pore_size.normal
pn.add_model(propname='pore.diameter',
             model=f,
             scale=1e-5, 
             loc=2.5e-5,
             seeds='pore.seed')
print(pn.models['pore.diameter@all'])
plt.hist(pn['pore.diameter'], edgecolor='k')

可以看到孔径分布并没有改成正态分布,这是因为pore.seed的值没有更改

print(pn.models['pore.seed@all'])
pn.models['pore.seed@all']['num_range'] = [0.01, 0.99]
pn.regenerate_models()
plt.hist(pn['pore.diameter'], edgecolor='k');

依赖项处理程序简介

孔隙尺度模型的主要优点是可以自动重新计算。例如,如果我们重新运行模型,那么依赖的所有其他模型也将自动重新计算

pn = op.network.Cubic(shape=[20, 20, 20], spacing=5e-5)
pn.add_model(propname='pore.seed',
             model=op.models.geometry.pore_seed.random)
pn.add_model(propname='pore.diameter',
             model=f,
             scale=1e-5, 
             loc=2.5e-5,
             seeds='pore.seed')
pn.add_model(propname='throat.seed',
             model=op.models.geometry.throat_seed.from_neighbor_pores,
             prop='pore.seed')
pn.add_model(propname='throat.diameter',
             model=f,
             scale=1e-5, 
             loc=2.5e-5,
             seeds='throat.seed')
print(pn['pore.seed'])
print(pn['throat.diameter'])
'''
[0.15319619 0.0328116  0.2287953  ... 0.23887453 0.8411705  0.10988224]
[6.59011212e-06 6.59011212e-06 1.75717988e-05 ... 1.13561876e-05
 2.66940036e-05 1.27284536e-05]
'''

pn.regenerate_models()
print(pn['pore.seed'])
print(pn['throat.diameter'])
'''
[0.36618201 0.35063831 0.40236746 ... 0.58822388 0.42461491 0.9016327 ]
[2.11640227e-05 2.11640227e-05 2.20171387e-05 ... 2.57748463e-05
 2.30989878e-05 3.79091080e-05]
'''

正如我们在上面看到的,调用前后的初始值都发生了变化。OpenPNM 通过检查每个函数的参数来自动执行此操作
下面的;流程图可以很好的解释依赖项处理程序的原理:

pn.models.dependency_map(style='planar');

绿色圆圈是孔隙属性,橙色方块是流道属性。我们还可以看到,有两个模型不依赖于其他任何东西:coordination_number和spacing。这些在生成过程中添加到网络模型,但不运行,因此它们的值在执行时不显示。

Everything not saved will be lost.
最后更新于 2023-03-30