物理应用基于matlab模拟井筒多相流【含Matlab源码 2152期】.zip

%=================================================================================================================================================================== % Sagd4LiHeel程序功能是:长管注汽(跟端定压力）水平段沿程参数计算 %=================================================================================================================================================================== function [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) [TubInfo]=Tubinfo(1); [Dtub,rto,htf,dertub]=deal(TubInfo(1,1),TubInfo(1,2),TubInfo(1,3),TubInfo(1,4)); % 长油管参数读入 [PWFN,TWFN,SQWFN,DPWF,DTWF,DSQWF,DPWFN,DTWFN,DSQWFN,TubPWFN,TubTWFN,TubSQWFN,TubLOSS,AunLOSS,QLOSS]=deal(zeros(1,NUM)); % 预分配 SIP=9.869; %% 传热方式及做功开关控制, 1 表示长油管到地层传热//摩擦力做功；2 表示长油管到环空传热//摩擦力不做功 [ft,R,Rtub,Ran]=Rselect(1,3,2,3,DELT); % ft及热阻计算选择 DQLctr=2; DWctr=2; %% 射孔摩阻控制, PerCtr=1,不算射孔摩阻;PerCtr=2,算射孔摩阻 PerCtr=2; % 算射孔摩阻如果算不了，可能是将割缝等效为射孔来计算射孔摩阻有问题，输出结果为复数，可尝试按射孔参数或者给摩阻系数赋定值来计算试试 %% 注汽井注入部分 %% 长油管沿程参数计算部分 QMctr=2; % 将地层温度设置高点，可以运行 QMctr=2 时的井筒沿程参数分布情况 if QMctr==1 QM=0; for M=1:NUM QM=QM+QG(M)*0.001; % 单位:吨/天 end elseif QMctr==2 QM=Vsx; % 水平井总流量,吨/天 end %% for M=1:NUM QMS=QM; % 长管注汽无质量损失，单位：吨/天 TT=T; DPWF(M)=-0.001; % 可设定不同的值对比结果差异，估算每个微元段的压力变化值，该值最好为负值，该处单位是atm DTWF(M)=-0.001; % 可设定不同的值对比结果差异，估算每个微元段的温度变化值，该值最好为负值，该处单位是℃ DSQWF(M)=-0.001; % 可设定不同的值对比结果差异，估算每个微元段的干度变化值，无因次 DPWFN(M)=inf; DTWFN(M)=inf; DSQWFN(M)=inf; while abs(DPWFN(M)-DPWF(M))>0.0001 || abs(DTWFN(M)-DTWF(M))>0.0001 || abs(DSQWFN(M)-DSQWF(M))>0.0001 if M==1 PWM=PWFC+DPWF(M)/2; % M微元段平均压力,单位是atm % PWFC(J)井跟端压力 TWM=TWFC+DTWF(M)/2; % M微元段平均压力,单位是℃ % TWFC(J)井跟端压力 SQWM=STQ+DSQWF(M)/2; % M微元段平均干度 % STQ(J)井跟端干度 elseif M~=1 PWM=TubPWFN(M-1)+DPWF(M)/2; % 单位是atm TWM=TubTWFN(M-1)+DTWF(M)/2; % 单位是℃ SQWM=TubSQWFN(M-1)+DSQWF(M)/2; end [~,~,~,~,~,~,DTP,ROSW,ROSM,HW,HS,DHST,DHWT,MIUW,MIUS]=STEAM(TWM,SQWM); [ROWM,ELL,Hl,vm,Rns,vsg]=Begtub(SQWM,QMS,ROSM,ROSW,Dtub); % 已核实 Reyn=vm*(Dtub/1000)*(ROSW*ELL+ROSM*(1-ELL))/(MIUW*ELL+MIUS*(1-ELL)); % 已核实，无因次，粘度单位是Pa.s,所以不需要跟产液井一样在vsl前面乘以1000 [f]=finj(5,Reyn,Dtub,dertub,ELL,Hl); twl=1/8*3.14159*f*Rns*(Dtub/1000)*DX*vm*vm; % 已核实 单位 N dwl=twl*vm; % 已核实 单位 功率 N*m/s=W=J/s tw=twl; dw=dwl; %% if DWctr==1 % 摩擦力做功 dw=dw+0; elseif DWctr==2 % 摩擦力不做功 dw=0; end %% PWM=PWM/SIP*1000000; % 单位是Pa DWL=dw/DX; % 已核实 单位为J/m.s %% if DQLctr==1 % 长油管到地层传热 DQL=(TWM-TT)/Rtub; % 单位为W/m,即J/m/s TubLOSS(M)=DQL*DX*86.4; % 已核实 转化为 KJ/d elseif DQLctr==2 % 长油管到环空传热,由于长油管温度大于环空，故环空向长油管传热不做考虑 DQL=(TWM-TWFN(M))/Ran; % 单位为W/m,即J/m/s end DPL=-f*Rns*vm^2/2/(Dtub/1000)/(1-ROWM*vm*vsg/PWM); % 已核实 DPL单位为Pa/m %% N1=QMS/86.4*(HS-HW)*1000; % 已核实 单位为J/s N2=QMS/86.4*(DHST*DTP*SIP/1000.0-DHWT*DTP*SIP/1000.0)*DPL; % 已核实 单位为J/m.s N3=DQL+DWL+QMS/86.4*DHWT*DTP*SIP/1000.0*DPL-QMS/86.4*vm*vsg/PWM*DPL; % 已核实 单位为J/m.s %% if M==1 DPWFN(M)=DPL*DX/1000000*SIP; % 单位转换为atm TubPWFN(M)=PWFC+DPWFN(M); % 每一微元段末尾端点处的压力值，单位atm TubSQWFN(M)=exp(-N2/N1*DX)*(-N3/N2*exp(N2/N1*DX)+STQ+N3/N2); % DX/2难道是井打在第一个网格中间其沿水平段压降第一网格只占1/2 if TubSQWFN(M)<=0.0 TubSQWFN(M)=0.0; end DSQWFN(M)=TubSQWFN(M)-STQ; end if M~=1 DPWFN(M)=DPL*DX/1000000*SIP; TubPWFN(M)=TubPWFN(M-1)+DPWFN(M); % PWFN(J,M)=PWFN(J,M-1)+0.5*(DPWFN(M)+DPWFN(M-1))；% 运算过程中SQWFN出现空值，原因是第一个循环假设的DPWF其实是0，这样就使得DPL为0，即分母N2为0 TubSQWFN(M)=exp(-N2/N1*DX)*(-N3/N2*exp(N2/N1*DX)+TubSQWFN(M-1)+N3/N2); if TubSQWFN(M)<=0.0 TubSQWFN(M)=0.0; end DSQWFN(M)=TubSQWFN(M)-TubSQWFN(M-1); end PPWF=TubPWFN(M)/SIP*1000; % 转换为千帕 if PPWF>=0.611 && PPWF<=22120 % 文献W.S.Tortlke,1989,Saturated Steam.Property Functional Correlations for Fully Implicit Thermal Reservoir Simulation T1=280.034+14.0856*log(PPWF)+1.38075*log(PPWF)*log(PPWF)-0.101806*log(PPWF)^3+0.019017*log(PPWF)^4; end TubTWFN(M)=T1-273.15; % 单位℃ if abs(DPWFN(M)-DPWF(M))<0.01 && abs(DSQWFN(M)-DSQWF(M))<0.01 break else DPWF(M)=DPWFN(M); DTWF(M)=DTWFN(M); DSQWF(M)=DSQWFN(M); end end %对应while end %对应for M=1:NUM %% 长油管注汽环空部分 for M=NUM:-1:1 QL=0.0; % 微元段水平井筒环空内流量 if M==NUM QL=0; elseif M~=NUM for N=M+1:NUM QL=QL+QG(N)*0.001; end end QMS=QM-QL; % 单位 吨/天 TT=T; DPWF(M)=-0.001; % 可设定不同的值对比结果差异，估算每个微元段的压力变化值，该值最好为负值，该处单位是atm DTWF(M)=-0.001; % 可设定不同的值对比结果差异，估算每个微元段的温度变化值，该值最好为负值，该处单位是℃ DSQWF(M)=-0.001; % 可设定不同的值对比结果差异，估算每个微元段的干度变化值，无因次 DPWFN(M)=inf; DTWFN(M)=inf; DSQWFN(M)=inf; while abs(DPWFN(M)-DPWF(M))>0.0001 || abs(DTWFN(M)-DTWF(M))>0.0001 || abs(DSQWFN(M)-DSQWF(M))>0.0001 if M==NUM PWM=TubPWFN(M)+DPWF(M)/2; % M微元段平均压力,单位是atm % PWFC(J)井跟端压力 TWM=TubTWFN(M)+DTWF(M)/2; % M微元段平均压力,单位是℃ % TWFC(J)井跟端压力 SQWM=T