二维微可压缩流体流动问题

油藏描述

本次算例计算一个二维油藏模型，网格数20*20，油藏大小100ft*100ft*5ft，四口井分布于（5，5）、（5，15）、（15，5）（15，15），四口井分别以100$ft^3/d$的产量生产。油藏边界为无流动边界，初始压力5500psi，油藏压缩系数为$1*10^{-7} 1/Pa$，初始孔隙度在0.3-0.4之间随机生成。x和y方向渗透率在0.6-0.4D之间随机生成。

流体描述

压力和粘度关系如下表所示

4000	5000	5500	6000	6500	7000	7500	8000
0.9100	0.9200	0.9243	0.9372	0.9650	0.9494	0.9812	1.0019

当计算压力值超出插值范围时使用表格中的最小值或者最大值进行替换。 初始油藏状态下的油体积系数为0.9，油压缩系数为$1*10^{-9} 1/Pa$。

方程求解

根据物质守恒方程、状态方程和达西定律，单相流动基本方程为

$$\frac{\partial }{\partial x}\left ( \beta_c \frac{A_x k_x}{\mu B} \frac{\partial p}{\partial x} \right )\Delta x + \frac{\partial }{\partial y}\left ( \beta_c \frac{A_y k_y}{\mu B} \frac{\partial p}{\partial y} \right )\Delta y - q_{sc }=\frac{V_b}{\alpha_c}\frac{\partial }{\partial t}\left ( \frac{\varphi }{B} \right )$$

$\beta_c$为传到率转换因子，在油田单位制下为1.127,在国际单位制下为1.

$\alpha_c$为体积转换系数，在油田单位制下为5.614583,在国际单位制下为1.

$V_b$为网格体积,$A$为截面积。

$$A_x = \Delta y * \Delta z$$ $$A_y = \Delta x * \Delta z$$

$\mu$是流体粘度，使用插值进行计算。

$B$为流体体积系数

$$B = \frac{B_0}{1 + c \left (p - p_0 \right )}$$

其中$B_0$是$p_0$时的体积系数，$c$为流体压缩系数.

$\varphi$为孔隙度

$$\varphi = \varphi_0 \left [1 + c_{\varphi } \left (p - p_0 \right ) \right ]$$

$c_{\varphi }$为油藏压缩系数.

基本方程变形为：

$$V_b \beta_c \frac{\partial }{\partial x}\left ( \beta_c \frac{k_x}{\mu B} \frac{\partial p}{\partial x} \right ) + V_b \beta_c \frac{\partial }{\partial y}\left ( \beta_c \frac{k_y}{\mu B} \frac{\partial p}{\partial y} \right ) - q_{sc }=\frac{V_b}{\alpha_c}\frac{\partial }{\partial t}\left ( \frac{\varphi }{B} \right )$$

使用有限差分逼近

$$T_{i + \frac{1}{2}, j} \left ( p_{i+1,j} - p_{i,j}\right ) - T_{i - \frac{1}{2}, j} \left ( p_{i,j} - p_{i-1,j}\right ) + T_{i, j + \frac{1}{2}} \left ( p_{i,j+1} - p_{i,j}\right )$$ $$- T_{i - \frac{1}{2}, j} \left ( p_{i,j} - p_{i-1,j}\right ) - q_{sc} = \Gamma_{i , j} \left( p_{i , j}^{n+1} - p_{i , j}^{n}\right )$$

其中

$$T_{i \pm \frac{1}{2},j} = \frac{2 T_{i,j} T_{i \pm 1,j} }{T_{i,j} + T_{i \pm 1,j}}$$ $$T_{i+1,j} = \frac{\Delta y \Delta z}{\Delta x} \beta_c \frac{k_x}{\mu B}_{i+1,j}$$ $$\Gamma_{i,j} = \left [ \frac{V_b}{\alpha_c \Delta t} \left ( \frac{\varphi c_\varphi}{B^n} + \frac{\varphi^n c}{B_0}\right )\right ]_{i,j}$$

由于$\frac{1}{\mu B}$是压力的函数，因此本方程为非线性方程。为了将方程线性化，使用上一个时间步的压力值计算T，其余地方的压力采用此时刻的值，因此本方法是采用半隐式方法进行求解。最终得到的半隐式方程为：

$$T_{i + \frac{1}{2}, j} p_{i+1,j} + T_{i - \frac{1}{2}, j} p_{i-1,j} + T_{i , j+ \frac{1}{2}} p_{i,j+1} + T_{i , j - \frac{1}{2}} p_{i,j-1} -$$ $$\left ( T_{i + \frac{1}{2}, j} + T_{i - \frac{1}{2}, j} + T_{i , j+ \frac{1}{2}} + T_{i , j - \frac{1}{2}} + \Gamma_{i,j}\right ) p_{i,j} = q_{sc} - \Gamma_{i,j} p_{i,j}^n$$

计算结果

采用时间步为60秒，开采20天后的压力分布结果如图

MatLab计算源码

%2D single phase with 4 production wells

%Basic Parameters, SI unit
clear all
dx = 5*0.3048;
dy = 5*0.3048;
dz = 5*0.3048;
nx = 20;
ny = 20;
beta_c = 1;
alpha_c = 1;
W = [5,5;5,15;15,5;15,15];      % well position
p0 = 5500*6894.76;      % original pressure 
por0 = 0.1*rand(nx*ny,1) + 0.3;
kx = (0.2*rand(nx*ny,1) + 0.4)*1e-12;       % random permeability in x direction
ky = (0.2*rand(nx*ny,1) + 0.4)*1e-12;       % random permeability in y direction
q = 200*0.3048^3/86400;     % m^3/s

p_mu = [4000;5000;5500;6000;7000;6500;7500;8000]*6894.76;
mu = [0.9100;0.9200;0.9243;0.9372;0.9650;0.9494;0.9812;1.0019]*1e-3;
B0 = 1.1;

c_phi = 1e-7;       % 地层压缩系数
c = 1e-9;           % 流体压缩系数

nt = 60*24*20;
dt = 60;

% number the grids
Grid = zeros();
N = 1;
for i = 1:nx
    for j = 1:ny
        Grid(j,i)= N;
        N = N+1;
    end
end

%Calculation
p = ones(nx*ny,nt)*p0;
T = zeros(nx*ny,nt);
G = zeros(nx*ny,nt);        % Gamma index
X = zeros(nx*ny,nx*ny);     % matrix
Y = zeros(nx*ny,1);        % right vector
u = ones(ny*nx,nt)*0.9E-3;
B = ones(ny*nx,nt)*1.1;
por = ones(ny*nx,nt)*0.3;
p_grid = zeros(ny,nx,nt);
p_grid(:,:,1) = reshape(p(:,1),ny,nx);

for t = 1:nt
    B(:,t) = B0./(1 + c_phi.*(p(:,t)-p0));   % fluid volume index
    por(:,t) = por0.*(1 + c_phi.*(p(:,t) - p0));
    T(:,t) = dz*beta_c*kx./u(:,t)./B(:,t);
    G(:,t) = dx*dy*dz/alpha_c/dt*(por0*c_phi./B(:,t) + por(:,t)*c./B0);
    
    Y = -G(:,t).*p(:,t);
    Y(Grid(W(1,1),W(1,2))) = q + Y(Grid(W(1,1),W(1,2)));
    Y(Grid(W(2,1),W(2,2))) = q + Y(Grid(W(2,1),W(2,2)));
    Y(Grid(W(3,1),W(3,2))) = q + Y(Grid(W(3,1),W(3,2)));
    Y(Grid(W(4,1),W(4,2))) = q + Y(Grid(W(4,1),W(4,2)));
    
    for i = 1:nx
        for j = 1:ny
            if p(Grid(i,j),t) >= 4000*6894.76 && p(Grid(i,j),t) <= 8000*6894.76
                u(Grid(i,j),t) = interp1(p_mu,mu,p(Grid(i,j),t),'linear')   ;  % viscosity
            elseif p(Grid(i,j),t) < 4000*6894.76
                u(Grid(i,j),t) = 0.91*1E-3;
            else
                u(Grid(i,j),t) = 1.0019*1E-3;
            end
            if mod(i,nx) == 1
                if mod(j,ny) == 1
                    X(Grid(j,i),((i-1)*ny+j)) = -2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t)) ...
                        -2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),i*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t));
                    X(Grid(j,i),((i-1)*ny+j+1)) = 2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) ;
                elseif mod(j,ny) == 0
                    X(Grid(j,i),((i-1)*ny+j)) =  -2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t))...
                        -2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),i*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t));
                    X(Grid(j,i),((i-1)*ny+j-1)) = 2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t));
                else
                    X(Grid(j,i),((i-1)*ny+j)) =  -2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t)) -...
                        2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t)) - ...
                        2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),i*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t));
                    X(Grid(j,i),((i-1)*ny+j-1)) = 2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t));
                    X(Grid(j,i),((i-1)*ny+j+1)) = 2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t));
                end
            elseif mod(i,nx) == 0
                if mod(j,ny) == 1
                    X(Grid(j,i),((i-1)*ny+j)) = -2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t)) ...
                        -2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),(i-2)*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t));
                    X(Grid(j,i),((i-1)*ny+j+1)) = 2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) ;
                elseif mod(j,ny) == 0
                    X(Grid(j,i),((i-1)*ny+j)) =  -2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t))...
                        -2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),(i-2)*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t));
                    X(Grid(j,i),((i-1)*ny+j-1)) = 2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t));
                else
                    X(Grid(j,i),((i-1)*ny+j)) =  -2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t)) -...
                        2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t)) - ...
                        2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) - G(Grid(j,i),t);
                   X(Grid(j,i),(i-2)*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t));
                    X(Grid(j,i),((i-1)*ny+j-1)) = 2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t));
                    X(Grid(j,i),((i-1)*ny+j+1)) = 2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t));
                end
            else
                if mod(j,ny) == 1
                    X(Grid(j,i),((i-1)*ny+j)) = -2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t)) - ...
                        2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t)) -...
                        2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),i*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t));
                    X(Grid(j,i),(i-2)*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t));
                    X(Grid(j,i),((i-1)*ny+j+1)) = 2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) ;
                elseif mod(j,ny) == 0
                    X(Grid(j,i),((i-1)*ny+j)) =  -2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t)) - ...
                        2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t)) -...
                        2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),i*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t));
                    X(Grid(j,i),(i-2)*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t));
                    X(Grid(j,i),((i-1)*ny+j-1)) = 2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t));
                else
                    X(Grid(j,i),((i-1)*ny+j)) =  -2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t)) - ...
                        2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t)) -...
                        2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t)) - ...
                        2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t)) - G(Grid(j,i),t);
                    X(Grid(j,i),i*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i+1),t)/(T(Grid(j,i),t)+T(Grid(j,i+1),t));
                    X(Grid(j,i),(i-2)*ny+j) = 2*T(Grid(j,i),t)*T(Grid(j,i-1),t)/(T(Grid(j,i),t)+T(Grid(j,i-1),t));
                    X(Grid(j,i),((i-1)*ny+j-1)) = 2*T(Grid(j,i),t)*T(Grid(j-1,i),t)/(T(Grid(j,i),t)+T(Grid(j-1,i),t));
                    X(Grid(j,i),((i-1)*ny+j+1)) = 2*T(Grid(j,i),t)*T(Grid(j+1,i),t)/(T(Grid(j,i),t)+T(Grid(j+1,i),t));
                end
            end
        end
    end
    
    p(:,t+1) = X\Y;
    p_grid(:,:,t+1) = reshape(p(:,t+1),ny,nx);
end

[m,n] = meshgrid(1:20);
surf(m,n,p_grid(:,:,end))       % final 3D plot