Orkiszewski算法的C语言实现

Orkiszewski的C语言实现中国石油大学（华东）石油工程学院Orkiszewski方法是计算气液两相垂直管流的方法，由于中间步骤涉及大量计算及迭代，使用计算机编程是实现计算的简便方法。笔者使用C语言编程。整个代码分为三大块，主程序、综合函数、功能函数。编程思路：多相管流计算所需的具体参数已经在程序中给出。采用从井底向上计算，按深度增量迭代。读者可以根据自己的需要修改参数和迭代方法。由于在一口井中一般很难同时出现四种流态。对于低含气井只会出现泡流或段塞流。下面所给程序中，没有给出过渡流的计算代码，请读者根据自身需要加以完善。编程中的主要难点是参数过多，所需计算的物性参数很多，导致定义变量时容易造成混乱。由于计算物性参数的程序和Orkiszewski程序的编写是分开进行的，讲两者很好的结合又是一个很大的问题。读者可以采用代码复用或文件包含的方法将二者结合。限于笔者的水平，以下所给的程序是将两者放在一个程序里，导致变量很多且混乱。请读者自行完善！如有好的建议，欢迎来信！参考资料：1.张琪.采油工程原理与设计.—东营：中国石油大学出版社2.步玉环、谷建伟、薛建泉.石油工程综合设计.—东营：中国石油大学出版社[e-mail]//奥奇斯泽斯基算法多相管流计算。质量含水率0.44，初始产液量设定为50t/d.//假设给定温度、压力、液体流速、流体密度、流体表面张力、管子直径、溶解气油比、生产气油比、压缩因子、原油体积系数//vs为滑脱速度：0.244m/s#includestdio.h#includemath.hdoubletempc(doubleQ,doubleH,doubleh,doubletr,doublet0,doubleFw);doubleRsc(doublegamao,doublegamag,doubleta,doublepa,doubleRp);doubleBoc(doubleRs,doublegamag,doublegamao,doubleta);doubledopc(doublegamao,doublegamag,doubleRs,doubleBo);doubleZc(doublegamag,doubleta,doublepa);doublevlc(doublegamaAPI,doubleta,doubleRs,doublefw);doublevgc(doublegamag,doublepa,doubleZ,doubleta);doublesigmac(doubleta,doublepa,doublegamaAPI,doublefw);doublezhhs(doublep1,intj);//定义综合函数。doubledz=100,H=2400;//定义两个全局变量，计算段长度和油层中深intmain(){charc;intj=0;doublep1;for(j=1;;j++){if(j==1){printf(Pleaseenterp1(MPa):);scanf(%lf,&p1);}p1=zhhs(p1,j);printf(计算点处深度：%.2lfm\n,H-j*dz);printf(Doyouwanttocontinue?(y/n):\n);getchar();//吸收回车c=getchar();if(c=='n'||c=='N')break;}printf(共进行了%d次，计算点处深度%lfm\n,j,H-j*dz);return0;}doublezhhs(doublep1,intj){/*q日产液量；Rp生产气油比；fw体积含水率；Fw质量含水率；dop油的密度；gamao油的相对密度；gamag气的相对密度；D管的直径（套管或油管）；tr油层温度；t3地表温度；dp设定压力增量；da大气密度；er误差精度；qo日产油量Qg日产气量；*/doubleq=50,Rp=20,fw=0.4,Fw,dop=850,gamao=0.85,gamag=0.7,D=0.124,tr=95,t3=20,dp=0.5,da=1.205,er=0.02;doubledl,dg,sigmal,sigma,gamaAPI,Rs,Z,ta,pa,p2,p2c=0,Bo,vl;doubleqo,Q,Qg,Wg,Wl,Ql,Qt,Wt,deta,vg,vs=0.244,vsg,vsi,Nb,taof,dpc;doubleAp,X,vt,vga,Hg,dm,vlh,Re,e,f;//D的值取决于是在套管中还是油管中doubleLB,LM,LS;intleap=0,lept=0,i=1;Fw=fw/(fw+(1-fw)*gamao);Ap=3.14/4*D*D;//产油量t/dqo=q*(1-Fw);loop:if(lept)dp=p1-p2c;p2=p1-dp;pa=(p2+p1)/2;pa*=1000000;//gamaAPIgamaAPI=141.5/gamao-131.5;//井筒温度场Q=q;ta=tempc(Q,H,j*dz,tr,t3,Fw);//溶解气油比Rs=Rsc(gamao,gamag,ta,pa,Rp);//体积系数Bo=Boc(Rs,gamag,gamao,ta);//原油密度dop=dopc(gamao,gamag,Rs,Bo);//混合液的密度dl=(dop+da*Rs*gamag)/Bo;printf(dl=%.2fkg/m3\n,dl);//压缩因子Z=Zc(gamag,ta,pa);//天然气密度dg=gamag*da*pa*293/Z/ta/0.101;//printf(dg=%.2fkg/m3\n,dg);//液体粘度vl=vlc(gamaAPI,ta,Rs,fw);//天然气粘度vg=vgc(gamag,pa,Z,ta);//气液表面张力sigma=sigmac(ta,pa,gamaAPI,fw);sigmal=sigma;//气体体积流量Qg=Z*qo*(Rp-Rs)*(ta+273.15)/(293.15*86400*gamao*pa*10);printf(Qg=%lfm3/s\n,Qg);//气体的质量流量Wg=1.206*qo*(Rp-Rs)*gamag/86400/gamao;printf(Wg=%lfkg/s\n,Wg);//液体的体积流量Ql=q*(Bo*(1-fw)+1*fw)/86400/(gamao*(1-fw)+1*fw);printf(Ql=%lfm3/s\n,Ql);//液体的质量流量Wl=(1000*q+1.206*gamag*Rs*qo/gamao)/86400;printf(Wl=%lfkg/s\n,Wl);//总体积流量Qt=Ql+Qg;printf(Qt=%lfm3/s\n,Qt);//总质量流量Wt=Wl+Wg;printf(Wt=%lfkg/s\n,Wt);//流型划分界限LB=1.071-0.7277*vt*vt/D;if(LB0.13)LB=0.13;vt=Qt/Ap;vga=Qg/Ap*pow(dl/9.81/sigmal,0.25);LS=50+36*vga*Ql/Qg;LM=75+84*pow(vga*Ql/Qg,0.75);if(Qg/QtLB)leap=1;//泡流if(Qg/QtLB&&vgaLS)leap=2;//段塞流if(vgaLS&&vgaLM)leap=3;//过渡流if(vgaLM)leap=4;//雾流if(leap==1)printf(泡流\n);elseif(leap==2)printf(段塞流\n);elseif(leap==3)printf(过渡流\n);elseif(leap=4)printf(雾流\n);//平均密度及摩擦损失梯度计算//泡流if(leap==1){//气相存容比Hg=0.5*(1+Qt/Ap/vs-sqrt(pow(1+Qt/Ap/vs,2)-4*Qg/vs/Ap));printf(Hg=%lf\n,Hg);//平均密度dm=dl*(1-Hg)+dg*Hg;printf(dm=%lf\n,dm);//液相真实流速vlh=Ql/Ap/(1-Hg);//雷诺数Re=dl*D*vlh/vl;//相对粗糙度e=2*4.57/100000/D;//摩擦阻力系数if(Re=2000)f=64/Re;elseif(Re=59.7/pow(e,8/7))f=0.3164/pow(Re,0.25);elseif(Re=(665-765*log10(e)/e))f=-1/1.8/log10(6.8/Re+pow(e/7.4,1.11));elsef=1/pow(2*log10(7.4/e),2);vlh=Ql/Ap/(1-Hg);taof=f*dl/D*vlh*vlh/2;}if(leap==2){//段塞流if(fw=0.5)//水是连续相{if(vt3.048)deta=0.00252*log10(1000*vl)/pow(D,1.38)-0.782+0.232*log10(vt)-0.428*log10(D);elsedeta=0.0174*log10(1000*vl)/pow(D,0.799)-1.352-0.162*log10(vt)-0.888*log10(D);}else{if(vt3.048)deta=0.00236*log10(1000*vl+1)/pow(D,1.415)-0.14+0.167*log10(vt)+0.113*log10(D);else{deta=0.00537*log10(1000*vl+1)/pow(D,1.371)+0.455+0.569*log10(D)-X;X=(log10(vt)+0.516)*(0.0016*log10(1000*vl+1)/pow(D,1.571)+0.722+0.63*log10(D));}}Re=vt*D*dl/vl;vs=(0.546+8.47/1000000*Re)*sqrt(9.81*D);Nb=vs*D*dl/vl;if(Nb8000&&Nb3000){vsi=(0.251+8.74/1000000*Re)*sqrt(9.81*D);vs=0.5*(vsi+sqrt(vsi*vsi+11170*vl/dl/sqrt(D)));Nb=vs*D*dl/vl;}if(Nb=8000){vs=(0.35+8.74/1000000*Re)*sqrt(9.81*D);Nb=vs*D*dl/vl;}dm=(Wt+dl*vs*Ap)/(Qt+vs*Ap)+deta*dl;taof=f*dl*vt*vt/2/D*((Ql+vs*Ap)/(Qt+vs*Ap)+deta);}if(leap==3){//过渡流省略}if(leap==4){//雾流Hg=Qg/(Qg+Ql);vsg=Qg/Ap;dm=(1-Hg)*dl+Hg*dg;Re=dg*vsg*D/vga;//相对粗糙度e=2*4.57/100000/D;//摩擦阻力系数if(Re=2000)f=64/Re;elseif(Re=59.7/pow(e,8/7))f=0.3164/pow(Re,0.25);elseif(Re=(665-765*log10(e)/e))f=-1/1.8/log10(6.8/Re+pow(e/7.4,1.11));elsef=1/pow(2*log10(7.4/e),2);vlh=Ql/Ap/(1-Hg);taof=f*dl/D*vlh*vlh/2;}dpc=(dm*9.81+taof)/(1-Wt*Qg/Ap/Ap/pa)*dz;//计算该段末端的压力p2c=p1-dpc/1000000;//比较压力增量的假设值与计算值，er=0.02;if(fabs(p2c-p2)er){i++;lept=1;gotoloop;}elseprintf(p2=%lfMPa,p2c=%lfMPa,dpc=%lfMPa