基于有限差分法-嵌入式离散裂缝网络的油气藏地层压力场计算：Matlab实现与理论实践指导

原创于 2026-03-13 10:15:00 发布 · 75 阅读

2 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

基于有限差分-嵌入式离散裂缝网络（FDM-EDFM）的油气藏地层压力场计算，通过matlab代码实现，可提供理论指导和相关问题，可计算不同裂缝网络的压力分布。

裂缝性油气藏的压力分布计算是个磨人的活儿。传统方法处理复杂裂缝网络就像用筷子夹豆腐——容易碎。今天咱们玩点实在的，用MATLAB把有限差分和嵌入式离散裂缝网络（FDM-EDFM）揉在一起，直接撸出压力场。

先看基质系统的处理。咱用结构化网格划分地层，这段代码生成网格坐标系：

Lx = 1000; Ly = 500; % 地层尺寸（米）
Nx = 50; Ny = 25;    % 网格数
dx = Lx/Nx; dy = Ly/Ny;
x = linspace(dx/2, Lx-dx/2, Nx);
y = linspace(dy/2, Ly-dy/2, Ny);
[xx, yy] = meshgrid(x, y);

这里用等距网格主要是为了方便后续处理裂缝的交切关系。注意网格中心点坐标的生成方式，这种半格距偏移能避免边界处理时踩坑。

裂缝处理是重头戏。先定义裂缝对象：

fractures = struct(...
    'points', {[200,100; 800,400], [500,50; 500,450]},... % 裂缝端点坐标
    'aperture', [0.1, 0.2],... % 开度（米）
    'perm', [100, 200]... % 渗透率（毫达西）
);

这里用结构体存储裂缝数据比用多个数组更清晰。注意裂缝是用端点定义的线段，后续需要自动计算与网格的交叉关系。

关键来了——计算裂缝网格连接关系：

[conn_matrix, trans] = deal([]);
for f = 1:length(fractures)
    [cross_grids, intersections] = find_intersections(fractures(f), dx, dy);
    % 自定义函数计算裂缝穿过的网格和交点
    for k = 1:size(cross_grids,1)
        i = cross_grids(k,1); j = cross_grids(k,2);
        trans = [trans; calc_trans(fractures(f), intersections(k,:), dx, dy)];
        % 计算传导率并存储
        conn_matrix = [conn_matrix; [sub2ind([Ny,Nx],j,i), f, trans(end)]];
    end
end

这段代码实现了裂缝网络与基质的耦合。find_intersections函数需要处理线段与网格的几何关系，这里没展开写（涉及射线交叉判断算法）。重点注意sub2ind的使用——把二维网格索引转为一维编号，这是构建系数矩阵的常用技巧。

组装整体矩阵时，稀疏矩阵能救命：

A = sparse(Nx*Ny + Nf, Nx*Ny + Nf); % Nf是裂缝段数
b = zeros(Nx*Ny + Nf, 1);

% 填充基质部分
for i = 1:Nx*Ny
    A(i,i) = -2*(1/dx^2 + 1/dy^2); % 中心系数
    % 处理相邻网格...
end

% 填充裂缝耦合项
for k = 1:size(conn_matrix,1)
    m = conn_matrix(k,1); f = conn_matrix(k,2);
    A(m,m) = A(m,m) - conn_matrix(k,3);
    A(m,Nx*Ny+f) = conn_matrix(k,3);
    A(Nx*Ny+f,m) = conn_matrix(k,3);
    A(Nx*Ny+f,Nx*Ny+f) = A(Nx*Ny+f,Nx*Ny+f) - conn_matrix(k,3);
end

这种分块矩阵结构把基质和裂缝压力统一求解。注意裂缝变量接在基质变量之后，用索引偏移实现。稀疏矩阵存储节省了至少90%的内存，不信你试试用full矩阵跑100x100网格。

最后求解和可视化：

p = A\b; % 简单粗暴直接求解
pressure_matrix = reshape(p(1:Nx*Ny), Ny, Nx);
contourf(xx, yy, pressure_matrix, 20, 'LineColor','none');
colormap(jet); colorbar;

实际生产代码得用迭代法求解，但教学演示直接反斜杠最直观。注意裂缝压力也可以单独提取分析：p(Nx*Ny+1:end)就是各裂缝段的压力值。

几个避坑指南：