一、源代码

function [Psum,h,V_gas,detap,P]=boreholepressure(Q_liquid)
Q_liquid=1;
%%%------------------------基础参数-----------------------------
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=1800  %%开始造斜
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
%%---------------------------计算井筒压力剖面--------------------------------------
 
 

三、备注

完整代码或者代写添加QQ912100926
往期回顾>>>>>>
【物理应用】杨氏双缝干涉实验【Matlab 310期】
【物理应用】井筒多相流【Matlab 311期】
【物理应用】二维对流扩散温度场【Matlab 312期】
【物理应用】内弹道【Matlab 313期】
【物理应用】两个光纤光栅的新光子微波发生器的仿真分析【Matlab 314期】