# [OpenPNM] Part 5——Examples of Algorithm —— Simulating Single Phase Transport

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

OpenPNM V2简易教程可参考：https://zhuanlan.zhihu.com/p/488306648

# 模拟单相流

## 创建网络模型

import numpy as np
import openpnm as op
import matplotlib.pyplot as plt

op.visualization.set_mpl_style()

pn = op.network.Demo(shape=[5, 5, 1], spacing=5e-5)
ax = op.visualization.plot_connections(pn)
ax = op.visualization.plot_coordinates(pn, ax=ax)
op.visualization.plot_tutorial(pn)

print(pn)


## 创建相

water = op.phase.Phase(network=pn)

#添加一个孔隙尺度模型来计算水的粘度：

model=op.models.phase.viscosity.water_correlation)
print(water)
print(water['pore.viscosity'])
'''
[0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
0.00089319 0.00089319 0.00089319 0.00089319 0.00089319 0.00089319
0.00089319]
'''


## 计算体积流率

### 简单的算法

Q=pi*R^4\Delta P/8\mu L

\Delta P 是压力损失

L 是细管长度

\mu 是黏度

Q 是体积流率

R 是半径

g_h=pi*R^4/8\mu L 称为水力传导系数，又称渗透系数

R = pn['throat.diameter']/2
L = pn['throat.length']
mu = water['throat.viscosity']  # See ProTip below
water['throat.hydraulic_conductance'] = np.pi*R**4/(8*mu*L)
print(water['throat.hydraulic_conductance'])
'''
[5.49820616e-15 9.99180736e-15 1.76321626e-14 5.07941658e-14
1.64674238e-14 1.56232587e-14 1.96205099e-14 3.29420230e-14
2.68272625e-15 2.41058903e-15 1.29513705e-15 1.37408974e-15
1.32009256e-15 1.26713904e-15 4.42946704e-15 3.66379252e-14
2.80882220e-14 1.85490397e-15 2.03568403e-15 8.86119497e-14
6.09330649e-15 9.96513487e-15 1.60549434e-14 3.04048800e-14
8.60638597e-14 1.96692761e-14 2.64071141e-15 4.56641017e-15
1.48844528e-15 1.28969372e-14 7.23297734e-15 1.24406579e-15
3.35705479e-15 1.68077772e-15 1.15437363e-14 7.47824866e-15
1.52303130e-15 1.53995733e-15 9.79881973e-14 3.60120595e-14]
'''


water.interpolate_data('throat.viscosity')
water['throat.viscosity']
'''
array([0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319,
0.00089319, 0.00089319, 0.00089319, 0.00089319, 0.00089319])
'''


OpenPNM为此提供了这样一种快捷方式，例如，如果请求一个不存在的流道属性，它将尝试获取孔隙的对应值并对值进行线性插值以生成流道的值数组。

### 严谨的算法

1. 考虑每个半孔和流道的渗透率。
2. 流道长度应该考虑球形孔体和圆柱形流道之间的重叠部分来仔细计算。
3. 孔隙的净截面积是通过孔中心和孔与流道的交点之间的积分来计算的
print(pn)
pn['throat.hydraulic_size_factors']
'''
[[2.74055386e-16 1.79111358e-17 3.27301354e-16]
[3.63255347e-16 2.26972347e-17 3.23512572e-16]
[2.88639845e-16 1.66334187e-17 2.70663772e-16]
[2.08378814e-16 9.05573092e-18 1.79279846e-16]
[2.41047346e-16 1.01860637e-17 1.88532351e-16]
[1.88532351e-16 1.14942402e-17 3.06377616e-16]
[6.91599038e-16 7.94384341e-17 7.33276170e-16]
[2.87560873e-16 8.96086740e-18 1.54282093e-16]
[3.39009315e-16 1.70105711e-17 2.57155332e-16]
[2.57155332e-16 1.52114636e-17 2.68854437e-16]
[2.91059527e-16 2.18393246e-17 4.13171768e-16]
[1.06381806e-16 1.03160463e-18 3.35319054e-17]
[2.79970479e-17 4.67608999e-19 2.38078927e-17]
[2.38078927e-17 6.17285071e-19 7.89928756e-17]
[5.50899795e-16 4.86887382e-17 5.24692673e-16]
[5.24692673e-16 4.98014766e-17 5.69207454e-16]
[6.44601336e-16 5.76268801e-17 5.53255336e-16]
[2.96667276e-16 1.23677349e-17 2.05130578e-16]
[2.05130578e-16 1.22530283e-17 2.91432658e-16]
[1.12444070e-16 1.44145783e-18 4.46657991e-17]
[2.74055386e-16 1.72960370e-17 3.04721477e-16]
[2.59224373e-16 1.05277663e-17 1.88532351e-16]
[3.23512572e-16 2.49701058e-17 4.28471886e-16]
[2.70663772e-16 2.00843657e-17 4.07871557e-16]
[1.62848722e-16 6.95389550e-18 1.54282093e-16]
[3.67078481e-16 2.83494786e-17 4.23444187e-16]
[1.88532351e-16 9.64767948e-18 2.11038523e-16]
[4.01213836e-16 2.14344968e-17 2.91059527e-16]
[7.70450700e-16 8.91544104e-17 7.47360177e-16]
[5.93659807e-17 8.29515339e-19 3.35319054e-17]
[1.01002324e-16 1.11930428e-18 3.70875827e-17]
[5.79797727e-17 5.47276971e-19 2.38078927e-17]
[2.91059527e-16 2.07575126e-17 3.80441333e-16]
[5.96817367e-16 5.15350675e-17 5.24692673e-16]
[3.35319054e-17 1.00672933e-18 1.01549457e-16]
[3.70875827e-17 1.21523141e-18 1.18381191e-16]
[2.38078927e-17 6.08744683e-19 7.66784021e-17]
[3.06100627e-16 1.25785692e-17 2.05130578e-16]
[5.24692673e-16 4.72033285e-17 5.25821568e-16]
[1.21860727e-16 1.49864288e-18 4.46657991e-17]]
'''


Q=F_h \Delta P /\mu=g_h \Delta P

F_h=pi*R^4/8 L

Q=\Delta P /(\mu/F_h ,_i+\mu/F_h ,_j+\mu/F_h ,_k)

F_h = water['throat.hydraulic_size_factors']
water['throat.hydraulic_conductance'] = (mu * (1/F_h).sum(axis=1))**(-1)
print(water['throat.hydraulic_conductance'])
'''
[1.79031699e-14 2.24355192e-14 1.66408558e-14 9.26774109e-15
1.04025141e-14 1.17150032e-14 7.27093840e-14 9.21045492e-15
1.70601726e-14 1.52639504e-14 2.16784149e-14 1.11005046e-15
5.05167063e-16 6.68542806e-16 4.61498082e-14 4.71552966e-14
5.40551328e-14 1.25652687e-14 1.24510347e-14 1.54419931e-15
1.72915552e-14 1.07498356e-14 2.46208856e-14 2.00152770e-14
7.15723217e-15 2.77400028e-14 9.84728144e-15 2.12921344e-14
8.08217352e-14 8.94098524e-16 1.20349374e-15 5.93476871e-16
2.06414110e-14 4.87077477e-14 1.08383083e-15 1.30441885e-15
6.59442665e-16 1.27744646e-14 4.47964989e-14 1.60429524e-15]
'''

#孔隙尺度模型
model=op.models.physics.hydraulic_conductance.generic_hydraulic)
print(water['throat.hydraulic_conductance'])
'''
[1.79031699e-14 2.24355192e-14 1.66408558e-14 9.26774109e-15
1.04025141e-14 1.17150032e-14 7.27093840e-14 9.21045492e-15
1.70601726e-14 1.52639504e-14 2.16784149e-14 1.11005046e-15
5.05167063e-16 6.68542806e-16 4.61498082e-14 4.71552966e-14
5.40551328e-14 1.25652687e-14 1.24510347e-14 1.54419931e-15
1.72915552e-14 1.07498356e-14 2.46208856e-14 2.00152770e-14
7.15723217e-15 2.77400028e-14 9.84728144e-15 2.12921344e-14
8.08217352e-14 8.94098524e-16 1.20349374e-15 5.93476871e-16
2.06414110e-14 4.87077477e-14 1.08383083e-15 1.30441885e-15
6.59442665e-16 1.27744646e-14 4.47964989e-14 1.60429524e-15]
'''


## 创建算法对象

op.algorithms.StokesFlow

sf = op.algorithms.StokesFlow(network=pn, phase=water)
#Definition : StokesFlow(name='stokes_?', *, phase, network)
#A subclass of LinearTransport to simulate viscous flow.
print(sf)


## 指定边界条件

sf.set_value_BC(pores=pn.pores('left'), values=100_000)
sf.set_rate_BC(pores=pn.pores('right'), rates=1e-10)
print(sf)


soln = sf.run()#soln is a dict with the quantity as the key (i.e. 'pore.pressure') and the solution as the value (i.e an ndarray).
print(soln)
#None


ax = op.visualization.plot_connections(pn)
ax = op.visualization.plot_coordinates(pn, ax=ax, color_by=water['pore.pressure'])
print(sf)
print(sf['pore.pressure'][pn.pores('right')])
'''
[146047.14279219 144668.45513766 131818.50627857 123987.42334449
153492.37334527]
'''


Everything not saved will be lost.