渗流模型的实现与解读

什么是渗流：

流体渗流是指流体在孔隙介质中的流动。
渗流状态，是指系统中出现一个大的集群，能够将这些集群节点的和邻居节点的边界打通、渗透（只考虑上下左右四个方向的邻居，集群即团簇）。

模型详解

首先，我们考虑 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 asd[xiabiao[s][0]][xiabiao[s][1]] = 1a = 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  #处于次大团簇中的点赋值为5000a[a == b2] = 10000  #处于最大团簇中的点赋值为10000a[(a != 5000) & (a != 10000)] = 0 #其余的点均归为0c3 = Counter(b.flatten())  # 统计当前最大和次大团簇的大小并打印p11 = float('%.2f' % p1) #保留两位小数

5.绘图：（图像均为动态效果，展示部分截图）
这是矩阵变化的动态图画法

plt.figure(figsize=(10,10))#创建画布
plt.clf()  # 清除上一幅图ax = plt.imshow(a,cmap=cmap) #画出矩阵aplt.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()