油藏描述
本次算例计算一个二维油藏模型，网格数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 |
方程求解
根据物质守恒方程、状态方程和达西定律，单相流动基本方程为
$$\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
Comments
Post a Comment