【物理应用】井筒多相流matlab源码

最新推荐文章于 2024-01-20 10:06:21 发布

matlab科研社

最新推荐文章于 2024-01-20 10:06:21 发布

阅读量14

点赞数

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/Matlab_dashi/article/details/140115428

版权

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  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.

【物理应用】井筒多相流matlab源码_物理应用

matlab科研社

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
【物理应用】井筒多相流matlab源码

登录后复制 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Analytic model of SAGD %% LI PENG...
复制链接

扫一扫