什么是渗流:
流体渗流是指流体在孔隙介质中的流动。
渗流状态,是指系统中出现一个大的集群,能够将这些集群节点的和邻居节点的边界打通、渗透(只考虑上下左右四个方向的邻居,集群即团簇)。
模型详解
首先,我们考虑 L × L的格子,我们一个一个地遍历白色格子,以概率为p1,p0去给这些格子染色。我们到第一个格子,就抛一枚硬币,假设这枚硬币正反面不均匀,正面出现的概率为p1,反面出现的概率为1-p1。
如果硬币出现正面,那么我就把当前格子染成黑色;如果硬币出现反面,那么我就不用管,让格子保留原来的白色。
然后考虑所有的黑色格子,当相邻的格子都为黑色时,我们就将其染成同一颜色,即黑色格子山下左右为黑色时,这几个格子成为一个团簇,染成同一种颜色。
不难发现,当p1的概率不断增大时,即1出现的次数增加,则黑色格子会越来越多,所构成的团簇也会越来越大,最终合并为一个团簇。
实验流程(附代码)
1.将所需要用到的库导入到python环境中:
# -*-coding:utf-8-*-
import numpy as np #用于创建随机数
import matplotlib.pyplot as plt #绘图
import sys #用于更改系统迭代次数
from collections import Counter #计数器
import time # 用于记录程序运行时间
import matplotlib as mpl
2.创建原始矩阵
biao = np.random.choice(10000, 10000, replace=False)#创建一万个随机数,用于存放下标
nums = np.random.choice([0, 1], size=10000, p=[1, 0]) #也可以使用np.zeros,此处便于后期调试
sd = nums.reshape(100, 100)#创建100x100的全零矩阵
xiabiao = []
for i in range(100):
for j in range(100):
xiabiao.append([i, j])#存放矩阵的随机下标,共10000个
#该函数用于随机向矩阵a中投入1,共10000次,即1出现的概率为0~1
def juzhen(s):
global a
sd[xiabiao[s][0]][xiabiao[s][1]] = 1
a = sd.copy()
return a
3.遍历矩阵,进行团簇分类
使用递归函数进行遍历,含有周期性边界条件。
#### 含周期性条件!!!
# 该函数用于递归遍历矩阵,进行团簇分类
# 方法为将同一团簇赋为同一值,不在团簇中的点保持为1
def digui(i, j):
if a[i][j] == 1:
a[i][j] = se
# 向下遍历
# if中的语句用于判断边界,即到达矩阵边缘时又跳回起点,上下左右均以相同方法运行
if i == len(a) - 1:
digui(0, j)
elif a[i + 1][j] == 1:
digui(i + 1, j)
# 向上遍历
if i == 0:
digui(len(a) - 1, j)
elif a[i - 1][j] == 1:
digui(i - 1, j)
# 向右遍历
if j == len(a[i]) - 1:
digui(i, 0)
elif a[i][j + 1] == 1:
digui(i, j + 1)
# 向左遍历
if j == 0:
digui(i, len(a[i]) - 1)
elif a[i][j - 1] == 1:
digui(i, j - 1)
注意,这里如果每生成一次矩阵就遍历一遍的话,整个程序需要运行至少三分钟,所以我采用每投入50个点进行一次团簇分类,这样就能节省很多时间:
if dd%50==0:
se = 2 #每个团簇有不同的颜色,即se的值不同
for i in range(len(a) - 1):
for j in range(len(a[i] - 1)):
if a[i][j] == 1:
digui(i, j)
se += 1 #每遍历完一个团簇就将颜色的值加1
4.统计最大和次大团簇
b = a.reshape(-1) #将矩阵转为一维数组,方便统计
b3 = b[b != 0] #去除0元素
b1 = np.bincount(b3) #统计数组中的团簇
b2 = np.argmax(b1) #取出最大团簇
b4 = b3[b3 != b2] #去除最大团簇
if b4.size != 0:
b11 = np.bincount(b4) #再次统计
b12 = np.argmax(b11) #取出次大团簇
a[a == b12] = 5000 #处于次大团簇中的点赋值为5000
a[a == b2] = 10000 #处于最大团簇中的点赋值为10000
a[(a != 5000) & (a != 10000)] = 0 #其余的点均归为0
c3 = Counter(b.flatten()) # 统计当前最大和次大团簇的大小并打印
p11 = float('%.2f' % p1) #保留两位小数
5.绘图:(图像均为动态效果,展示部分截图)
这是矩阵变化的动态图画法
plt.figure(figsize=(10,10))#创建画布
plt.clf() # 清除上一幅图
ax = plt.imshow(a,cmap=cmap) #画出矩阵a
plt.title("P1:%f"%p11)
plt.xlabel("Maximum cluster(red):%d "
"Secondary cluster(blue):%d"%(c3[10000],c3[5000]))
plt.pause(0.1) # 暂停一段时间,不然画的太快会卡住显示不出来
这是最大团簇的变化趋势的折线图画法:
图中柱状图为最大团簇的大小,折线图为最大团簇的变化率
plt.close()
fig = plt.figure('size & change rate',figsize=(10, 6)) # 整体图的标题
plt.bar(p111, zui,[0.002],color='red')
plt.plot(p111, bian)
plt.title("Maximum cluster&change rate")
plt.xlabel("Maximum cluster(bar) ""change rate")
plt.show()