boreholepressure_井底压力计算_井筒_井筒压力_

function [Psum,h,V_gas,detap,P]=boreholepressure(Q_liquid) %%%------------------------基础参数----------------------------- p_stand=0.101325; T_stand=273.15; %热力学温度 density_gas_stand=1.25/1000; V_stand=0.022414; T0=293; Gt=2.84; %%%%———————————————井身结构——————————————— H=43860; H_1=1850 %%开始造斜 H_2=2634 %造斜结束 d_bore =244.5; %%技术套管 d_piple=127; d_bore_1=215.9;%%三开井眼直径 d_collar=127; %177.8mm L_collar=150; %-----------------------钻井参数-------------------------------------- %Q_liquid=25; Q_gas=111; visocity_liquid=40; %%mpa s visocity_gas=25; Yp=8; %pa density_mud=1.48; deta_H=20; N=H/deta_H; Tg=zeros(N+1,1); %------------------注气量---------------- Ratio_gas=Q_gas/Q_liquid; p_gas=zeros(N+1,1);V_gas=zeros(N+1,1);;density_mix=zeros(N+1,1); visocity_mix=zeros(N+1,1);Q_mix=zeros(N+1,1);detap=zeros(N+1,1); v=zeros(N+1,1);Re=zeros(N+1,1);h=zeros(N,1); %----------------------生产井眼温度和静液柱压力剖面---------------------------------- for i =1:N+1 if i>1 && i<=H_1/deta_H+1 Tg(i,1)=T0+(i-1)*deta_H*Gt/100; p_gas(i,1)=(i-1)*deta_H*9.8*density_mud/1000; %%单位兆帕MPa elseif i>H_1/deta_H+1 && i<=H_2/deta_H+1 Tg(i,1)=Tg(H_1/deta_H+1,1)+39.2*(i-(H_1/deta_H+1))*Gt/100; p_gas(i,1)=9.8*density_mud/1000*H_1+12.7*(i-(H_1/deta_H+1))*9.8*density_mud/1000; %%单位兆帕MPa V_gas(i,1)=1/(1+p_gas(i,1)*T_stand/(Ratio_gas/p_stand/Tg(i,1))); elseif i>H_2/deta_H+1 && i<=N+1 Tg(i,1)=T0+(H_1+636)*Gt/100; p_gas(i,1)=(H_1+636)*9.8*density_mud/1000; %%单位兆帕MPa end end %---------------------计算气体浓度和混合物密度--------------------------- for i=1:N+1 V_gas(i,1)=1/(1+p_gas(i,1)*T_stand/(Ratio_gas*p_stand*Tg(i,1))); density_mix(i,1)=density_mud*(1-V_gas(i,1))+V_gas(i,1)*density_gas_stand*p_gas(i,1)*T_stand/Tg(i,1)/p_stand; end %---------------------计算混合物粘度及混合物粘度--------------------------------------- for i=1:N+1 if V_gas(i,1)>0 && V_gas(i,1)<=0.54 visocity_mix(i,1)=visocity_liquid*(1+3.6*V_gas(i,1)); elseif V_gas(i,1)>0.54 && V_gas(i,1)<=0.96 visocity_mix(i,1)=visocity_liquid/(1-V_gas(i,1)^0.49); else visocity_mix(i,1)=visocity_gas*(1+2.5*(1-V_gas(i,1))); end Q_mix(i,1)=(1+Ratio_gas*p_stand*Tg(i,1)/T_stand/p_gas(i,1))*Q_liquid; end %-------------------------计算环空返速与雷诺数-------------------------- for i=1:N+1 if i>0 && i<=H_2/deta_H+1 v(i,1)=Q_mix(i,1)/(d_bore^2-d_piple^2)/3.14*4*1000; %%%%单位 m/s; Re(i,1)=density_mix(i,1)*v(i,1)*(d_bore-d_piple)/visocity_mix(i,1)*1000; elseif i>H_2/deta_H+1 && i<=(H-L_collar)/deta_H+1 v(i,1)=Q_mix(i,1)/(d_bore_1^2-d_piple^2)/3.14*4*1000; Re(i,1)=density_mix(i,1)*v(i,1)*(d_bore_1-d_piple)/visocity_mix(i,1)*1000; else v(i,1)=Q_mix(i,1)/(d_bore_1^2-d_collar^2)/3.14*4*1000; Re(i,1)=density_mix(i,1)*v(i,1)*(d_bore_1-d_collar)/visocity_mix(i,1)*1000; end end %---------------------------计算摩阻系数------------------------------ f=zeros(N+1,1); for i=1:N+1 if i>0 && i<=H_2/deta_H+1 ypsl=2*0.5/d_bore; if Re(i,1)<=2000 %%层流区 f(i,1)=64/Re(i,1); elseif Re(i,1)>2000 && Re(i,1)<=59.7/ypsl^(8/7) %%水力光滑区 f(i,1)=0.3164/(Re(i,1)^0.25); elseif Re(i,1)>59.7/ypsl^(8/7) && Re(i,1)<=(665-765 *log(ypsl))/ypsl %%混合摩擦区 f(i,1)=1/(-1.8*(log(6.8/Re(i,1)+(0.5/(3.7 * d_bore))^1.11)/log(10)))^2; elseif Re(i,1)>((665-765*(log(ypsl)/log(10)))/ypsl) f(i,1)=1/(2*(log(3.7*d_bore/0.5)/log(10)))^2; end else ypsl=2*1/d_bore_1; if Re(i,1)<=2000 %%层流区 f(i,1)=64/Re(i,1); elseif Re(i,1)>2000 && Re(i,1)<=59.7/ypsl^(8/7) %%水力光滑区 f(i,1)=0.3164/(Re(i,1)^0.25); elseif Re(i,1)>59.7/ypsl^(8/7) && Re(i,1)<=(665-765 *log(ypsl))/ypsl %%混合摩擦区 f(i,1)=1/(-1.8*(log(6.8/Re(i,1)+(1/(3.7 * d_bore_1))^1.11)/log(10)))^2; elseif Re(i,1)>((665-765*(log(ypsl)/log(10)))/ypsl) f(i,1)=1/(2*(log(3.7*d_bore_1/1)/log(10)))^2; end end end %------------------------计算环空摩擦压降---------------------- for i=1:N+1 if i>0 && i<=H_2/deta_H+1 detap(i,1)=2*f(i,1)*density_mix(i,1)*deta_H*v(i,1)^2/(d_bore-d_piple); elseif i>H_2/deta_H+1 && i<=(H-L_collar)/deta_H+1 detap(i,1)=2*f(i,1)*density_mix(i,1)*deta_H*v(i,1)^2/(d_bore_1-d_piple); else detap(i,1)=2*f(i,1)*density_mix(i,1)*deta_H*v(i,1)^2/(d_bore_1-d_collar); end end %%---------------------------计算井筒压力剖面-------------------------------------- detaP=detap(2:N+1);Density_mix=density_mix(2:N+1);Psum=zeros(N,1); for i =1:N h(i,1)=i*deta_H; end for i=1:N if i>0 && i<=H_1/deta_H+1 Psum(i,1)=sum(detaP(1:i))+0.00981*Density_mix(i,1)*deta_H*i+4; elseif i> H_1/deta_H+1 && i<=H_2/deta_H+1 Psum(i,1)=sum(detaP(1:i))+0.00981*Density_mix(i,1)*H_1+12.7*0.00981*Density_mix(i,1)*(i-H_1/deta_H-1)+4; else Psum(i,1)=sum(detaP(1:i))+0.00981*Density_mix(i,1)*H_1+0.00981*Density_mix(i,1)*636+4; end end P=zeros(N,1); for i=2:N P(i,1)=sum(detap(2:i)); end