%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analytic model of SAGD %
% LI PENG %
% China university of petroleum %
% DATE 2019/5/4 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 初始化Matlab,工作空间清理
clc;
clear;
tic;
%% 数据载入
SIP=9.869; % 读入参数 % SIP:大气压atm与MPa之间的换算
[Info]=Info(1); % 水平井水平段信息
[WellSlotScr]=WellSlot(1);
%% 参数读取
[DX,DY,DZ]=deal(Info(1,1),Info(1,2),Info(1,3)); % (1)网格数据读入
[P,T,SO,SW,KX,KY,KZ,VP]=deal(Info(2,1),Info(2,2),Info(2,3),Info(2,4),Info(2,5),Info(2,6),Info(2,7),Info(2,8)); %(2)油藏初始数据读入
[PWFC,TWFC,STQ,S,L,Vsx,PWFCS,TWFCS,STQS]=deal(Info(3,1),Info(3,2),Info(3,3),Info(3,4),Info(3,5),Info(3,6),Info(3,7),Info(3,8),Info(3,9)); %(3)井控数据读入
[D,rlo,rh,hf,ha,lamdace,lamdae,arfa,ls,ws,lu,ts,reh,dert,ms,ns,ngf]=deal(WellSlotScr(1,1),WellSlotScr(1,2),... %(4)筛管数据读入
WellSlotScr(1,3),WellSlotScr(1,4),WellSlotScr(1,5),WellSlotScr(1,6),WellSlotScr(1,7) ,WellSlotScr(1,8),...
WellSlotScr(1,9), WellSlotScr(1,10),WellSlotScr(1,11),WellSlotScr(1,12),WellSlotScr(1,13),WellSlotScr(1,14),...
WellSlotScr(1,15),WellSlotScr(1,16),WellSlotScr(1,17));
%% 统一压力单位
P=P*SIP; %※压力相关的要单位转为 atm
PWFC=PWFC*SIP; %※压力相关的要单位转为 atm
PWFCS=PWFCS*SIP; %※压力相关的要单位转为 atm
SG=1-SO-SW;
%%
NUM=L/DX; % 水平段微元段数目
DELT=1; % 时间步长,day
FT=0; % 已模拟总时间,day
FTMAX=1.5; % 模拟总时间,day
%%
TOC=zeros(1,1); % 总产油量
TWC=zeros(1,1); % 总注/产水量
Matinj1=[ ];
Matinj2=[ ];
Matinj3=[ ];
%% 大循环
NNMAX=8000000;
Iteration=0; % 迭代次数记录
for N=1:NNMAX
FT=FT+DELT;
if FT<=FTMAX
PWF(1:NUM)=PWFC;
TWF(1:NUM)=TWFC;
SQWF(1:NUM)=STQ;
elseif FT>FTMAX
break
end
[DDP,DDSQ]=deal(1,0);
while DDP>=0.01 || DDSQ>0.01 % while判断循环:井筒压力、干度
Iteration=Iteration+1; % 迭代次数记录
[QO,QW,QG,QHG]=QRATEsteam(NUM,P,T,SW,SG,S,PWF,TWF,SQWF,DX,DY,DZ,KY,KZ,rlo); %※ 注汽解析模型
%% 注汽解析模型部分
[PWFN,TWFN,SQWFN,TubPWFN,TubTWFN,TubSQWFN,TubLOSS,AunLOSS,QLOSS]=Sagd4LiHeel(NUM,DX,Vsx,QG,T,PWFC,TWFC,STQ,DELT,D,ls,ws,lu,dert,ngf); % 可运算
%% 注汽判断
DDP=max(abs((PWF(1:NUM)-PWFN(1:NUM))),[],2); % 求出每行的最大值
DDT=max(abs((TWF(1:NUM)-TWFN(1:NUM))),[],2); % PWFN和TWFN为wellbore中计算最新值;应该是计算值与假定值的对比
DDSQ=max(abs((SQWF(1:NUM)-SQWFN(1:NUM))),[],2);
PWF(1:NUM)=PWFN(1:NUM);
TWF(1:NUM)=TWFN(1:NUM);
SQWF(1:NUM)=SQWFN(1:NUM);
if DDP<10 && DDSQ<10 % 该处判断其实可以去掉。误差值0.0001是否太小,调大点可以减少循环次数
break;
end
%% 迭代次数达到一定值跳出循环,注意迭代次数越大越接近精确值
Iteration
if Iteration>=1000
break
end
%%
end % 对应while DDP>=0.0001 || DDT>=0.0001 || DDSQ>0.0001
%% 井筒压力单位处理
TubPWFN(:)=TubPWFN(:)/SIP; % 输出单位为MPa
PWFN(:)=PWFN(:)/SIP; % 输出单位为MPa
%% 注汽解析模型数据输出
TWR=0;
for M=1:NUM
PP=P;
TT=T;
TIN=TWFN(M);
TIT=(TT+TIN)/2.0; % ???是否加权平均
[RG,~,HG]=STEAM(TWFN(M),SQWFN(M));
[TubRG,~,TubHG]=STEAM(TubTWFN(M),TubSQWFN(M));
HG=HG/1000; % 焓值单位 10^6 J/kg
TubHG=TubHG/1000; % 焓值单位 10^6 J/kg
QGG=QG(M)/1000; % QG单位为Kg/D,该处除以水密度则为 m3/D,按地面条件下水的体积来算
QtonD=QG(M)/1000; % 注汽量输出单位为 吨/天
QKGS=QG(M)/86400; % 注汽量输出单位为 kg/s
Entha=QHG(M)/DX/86400; % QHG单位是 kJ/m.s
Tubloss=TubLOSS(M)/DX/86400; % 单位是 kJ/m.s
Aunloss=AunLOSS(M)/DX/86400; % 单位是 kJ/m.s
Qloss=QLOSS(M)/DX/86400; % 单位是 kJ/m.s
HeatTA=Entha+Tubloss+Aunloss; % 单位是 kJ/m.s ;注入热量加长油管/短油管到地层以及环空到地层的热量;根据井筒中传热方式计算选择地层到底吸收了多少热量
HeatA=Entha+Aunloss; % 单位是 kJ/m.s ;注入热量加环空到地层的热量;根据井筒中传热方式计算选择地层到底吸收了多少热量
Matinj1=[Matinj1;FT,M,M*DX,TubPWFN(M),PWFN(M),TubTWFN(M),TWFN(M),TubSQWFN(M),SQWFN(M),QGG,QtonD,QKGS]; % Matinj1指输出的第1个注入井数据矩阵
Matinj2=[Matinj2;FT,M,M*DX,TubHG,HG,Tubloss,Aunloss,Qloss,Entha,HeatTA,HeatA]; % 之前输出顺序为[TubHG,HG,Tubloss,Aunloss,HeatTA,HeatA,Entha,Qloss]
TWR=TWR+QtonD; % 单位为 吨/天
TWC=TWC+QtonD*DELT; % 单位为 吨
end
Matinj3=[Matinj3;FT,DELT,TWR,TWC]; % Matpro2指输出的第2个注入井数据矩阵
%%
N
toc
end % 对应 for N=1:NNMAX
%% 出图部分
% hold on
% figure;plot(Matinj1(:,3),Matinj1(:,[4,5]),'r+') figure;plot(Matinj1(:,3),Matinj1(:,4),'r+')
% figure;plot(Matinj1(:,3),Matinj1(:,4),'LineStyle', 'none', 'Marker', '^')
figure;plot(Matinj1(:,3),Matinj1(:,4),'r^')
hold on
plot(Matinj1(:,3),Matinj1(:,5),'bs')
title('井筒蒸汽压力分布图','Fontsize',18)
xlabel('距水平井跟端距离/m','Fontsize',17)
ylabel('井筒蒸汽压力/MPa','Fontsize',17)
h=legend('长管','环空');
set(h,'FontSize',13); %
set(gca,'FontSize',11); % 坐标字体大小设置
%
figure;plot(Matinj1(:,3),Matinj1(:,6),'r^')
hold on
plot(Matinj1(:,3),Matinj1(:,7),'bs')
title('井筒蒸汽温度分布图','Fontsize',18)
xlabel('距水平井跟端距离/m','Fontsize',17)
ylabel('井筒蒸汽温度/°C','Fontsize',17)
h=legend('长管','环空');
set(h,'FontSize',13); %
set(gca,'FontSize',11); % 坐标字体大小设置
%
figure;plot(Matinj1(:,3),Matinj1(:,8),'r^')
hold on
plot(Matinj1(:,3),Matinj1(:,9),'bs')
title('井筒蒸汽干度分布图','Fontsize',18)
xlabel('距水平井跟端距离/m','Fontsize',17)
ylabel('井筒蒸汽干度/°C','Fontsize',17)
h=legend('长管','环空');
set(h,'FontSize',13); %
set(gca,'FontSize',11); % 坐标字体大小设置
%% 数据输出
% Nameinj1={'FT','M','TubPWFN(M)','PWF(M)','TubTWFN(M)','TWF(M)','TubSQWFN(M)','SQWF(M)','QGG','QKGS','RO','RW'}; % 注入井数据输出
% xlswrite('Parameters',Nameinj1,'sheet1','A1');
% xlswrite('Parameters',Matinj1,'sheet1','A2');
%
% Nameinj2={'FT','DELT','日注水TWR','总注水TWC(J)'};
% xlswrite('Parameters',Nameinj2,'sheet1','R1');
% xlswrite('Parameters',Matinj2,'sheet1','R2');
toc
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
- 12.
- 13.
- 14.
- 15.
- 16.
- 17.
- 18.
- 19.
- 20.
- 21.
- 22.
- 23.
- 24.
- 25.
- 26.
- 27.
- 28.
- 29.
- 30.
- 31.
- 32.
- 33.
- 34.
- 35.
- 36.
- 37.
- 38.
- 39.
- 40.
- 41.
- 42.
- 43.
- 44.
- 45.
- 46.
- 47.
- 48.
- 49.
- 50.
- 51.
- 52.
- 53.
- 54.
- 55.
- 56.
- 57.
- 58.
- 59.
- 60.
- 61.
- 62.
- 63.
- 64.
- 65.
- 66.
- 67.
- 68.
- 69.
- 70.
- 71.
- 72.
- 73.
- 74.
- 75.
- 76.
- 77.
- 78.
- 79.
- 80.
- 81.
- 82.
- 83.
- 84.
- 85.
- 86.
- 87.
- 88.
- 89.
- 90.
- 91.
- 92.
- 93.
- 94.
- 95.
- 96.
- 97.
- 98.
- 99.
- 100.
- 101.
- 102.
- 103.
- 104.
- 105.
- 106.
- 107.
- 108.
- 109.
- 110.
- 111.
- 112.
- 113.
- 114.
- 115.
- 116.
- 117.
- 118.
- 119.
- 120.
- 121.
- 122.
- 123.
- 124.
- 125.
- 126.
- 127.
- 128.
- 129.
- 130.
- 131.
- 132.
- 133.
- 134.
- 135.
- 136.
- 137.
- 138.
- 139.
- 140.
- 141.
- 142.
- 143.
- 144.
- 145.
- 146.
- 147.
- 148.
- 149.
- 150.
- 151.
- 152.
- 153.
- 154.
- 155.
- 156.
- 157.
- 158.
- 159.
- 160.
- 161.
- 162.
- 163.
- 164.
- 165.
- 166.
- 167.
- 168.
- 169.
- 170.
- 171.
- 172.
- 173.
- 174.
- 175.
- 176.
- 177.
- 178.
- 179.
- 180.
- 181.
- 182.