Matnote_1_MATLAB实现渗透(渗滤)模型 Percolation Model

hyjack-00

于 2022-01-26 00:12:54 发布

阅读量2.3k

点赞数 9

分类专栏： Matnote 文章标签： matlab 开发语言算法

本文链接：https://blog.csdn.net/hyjiesix/article/details/122693406

版权

Matnote 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

在MATLAB中采用蒙特卡罗方法仿真模拟并可视化物理渗透模型，具体采用快速并查集算法实现动态连通系统，并利用平均思想估计渗透阈值常数。

问题

实际物理问题：已知一多孔材料的孔隙率，求问液体能否从顶端贯穿该材料直至到达底部；有一个由绝缘和金属材料组成的复合系统，求问需要有多少部分材料是金属的才能使系统成为电导体……

我们将其此类宏观连通性问题抽象为数学渗透模型：在 N*N 个站点(site) 的网格中每个站点有开(open) 与闭(blocked) 两种状态，若存在一条开路连通网格顶层与底层，则称该网格渗透(percolate)。我们要研究的是网格的渗透概率与站点空缺率p 之间的关系。

请添加图片描述
显然随着空缺站点占比的增加，网格渗透的可能性随之增加。研究显示随着网格尺寸 N 的增加，这种变化趋势越来越趋向于一种相变。在 N 充分大时，有一个常数 P 使得 p < P 时随机网格几乎从不渗透，而 p > P 时随机网格总是渗透。

请添加图片描述

这个常数 P 称为渗透阈值(Percolation threshold)，其值约为0.592746。目前尚无确定渗透阈值的数学解决方案，我们的任务是通过计算机仿真模拟来复现这个现象并估算渗透阈值 (￣▽￣)。

思路

蒙特卡罗：重复以下步骤：
- 初始化一个全部闭锁的网格
- 每次随机增加一个开站点
- 判断顶层与底层是否连通(渗透)
- 如果连通，则记录下此时的站点空缺率，进行下一轮初始化；如果不连通则继续开启站点
动态连通：使用路径压缩的加权 Union-Path 算法，每一次开启新站点，将其与上下左右的开站点依次进行并(union) 操作。判断是否渗透即询问顶层与底层是否连通，这里可以增加 top 和 bottom 两个虚拟站点来进行询问。
渗透阈值：记录连通瞬间的空缺率，绘出关于空缺率的渗透频数图像，其分布应该会大致类似于正态分布，再关于空缺率对频数进行归一化与累加就应当得到上一张图中的相变曲线。估算的渗透阈值可以直接由频数分布的平均值得到。

模拟结果

MATLAB画图是真不戳~

这是不同测试条件下的模拟结果，从左上到右下分别有：网格尺寸 N 取 10, 20, 50, 100，测试次数取 10k, 1k, 1k, 1k。
请添加图片描述
对其进行累加操作得到越来越趋向于相变的曲线，验证成功：

图中绿线处即为渗透阈值。四种情况下阈值的平均值分别为 0.5900, 0.5905, 0.5921, 0.5912，与实际的阈值常数有些许误差，时间有限不在这里分析误差了ヽ(ー_ー)ノ。

可视化

不要在意渗透的方向（画出来才知道是横着的。。）

N = 20 （动图为两次模拟）

在这里插入图片描述
N = 100 （动图为一次模拟）

代码

数据处理与可视化的部分省略

%% Percolation 
% Basic parameter
    
global N Grid id sz 
test  = 3;  % 30 times: big enough
N     = 20;
pcl  = zeros(1001, 1);          % Percolation Frequency
px   = linspace(0, 100, 1001);  % Site Vacancy Rate \%

%% 1:blocked  2:open(connected to top)   3:open
% name site id from 1 to N*N+2 (1:top  N*N+2:bottom)
fms  = 1;


for freq = 1:test
    G    = zeros(N, N);
    Grid = ones(N*N, 1);
    open = 0;
    id   = ones(N*N+2, 1);     % id pointer
    sz   = ones(N*N+2, 1);     % size of tree
    for i = N+1:N*N+1
        if i > N*N-N
            id(i+1) = N*N+2;
        else
            id(i+1) = i+1;
        end
    end
    
    while 1 % Monte Carlo
        x = fix(N * N * rand) + 1;   % random blocked site x
        if Grid(x) ~= 1
            continue
        end
        Grid(x) = 3;
        open = open + 1;

        x0 = fix((x-1)/N)+1;
        y0 = x-x0*N+N;
        root1 = find(x+1);   % blocked site
        if y0 > 1   % left
            if Grid(x-1) ~= 1
                root2 = find(x);
                root1 = union(root1, root2);
            end
        end
        if x0 > 1   % up
            if Grid(x-N) ~= 1
                root2 = find(x-N+1);
                root1 = union(root1, root2);
            end
        end
        if y0 < N   % right
            if Grid(x+1) ~= 1
                root2 = find(x+2);
                root1 = union(root1, root2);
            end
        end
        if x0 < N   % down
            if Grid(x+N) ~= 1
                root2 = find(x+N+1);
                root1 = union(root1, root2);
            end
        end

        if  find(1) == find(N*N+2)   % percolation judging
            temp = fix(open*1000/(N*N));
            pcl(temp) = pcl(temp) + 1;
            break
        end
    end
end

function Uroot = union(root1, root2)
    global id sz
    if root1 == root2
        Uroot = root2;
    elseif sz(root1) > sz(root2)
        sz(root1) = sz(root1) + sz(root2);
        id(root2) = root1;
        Uroot = root1;
    else
        sz(root2) = sz(root1) + sz(root2);
        id(root1) = root2;
        Uroot = root2;
    end
    % ignore the same root
end

function root = find(p)
    global id
    while p ~= id(p)
        id(p) = id(id(p));
        p = id(p);
    end
    root = p;
end