#include "math.h"
#include "stdio.h"
#define  GOR  15       /****氣油比****/
#define  m    0.033     /****地溫梯度****/
#define  t0   90      /****地層原始溫度****/
#define  y    24        /****日產(chǎn)油量****/
#define  s    15        /****工作周期****/
#define  cp1  2.082     /****原油比熱****/
#define  r1   0.9727    /****20℃下原油相對密度****/
#define  cp2  2.4       /****伴生氣比熱****/
#define  r2   0.6154    /****伴生氣密度****/
#define  xs   0.81216   /****水泥導熱系數(shù)****/
#define  xd   1.74528   /**** 地層導熱系數(shù)****/
#define  a    1.032e-6  /****地層熱擴散系數(shù)****/
#define  C    33.85    /****高凝油的含蠟量****/
#define  qi   55        /****電加熱功率****/
#define  PI   3.1415    /****圓周率****/
#define  g    9.8       /****重力加速度****/
 double niandu(double t)  /******t溫度地面脫氣原油的粘度******/
 { 
	 double x,nd;
   x=-4.009*log10(t)+10.0942;
   nd=exp(x)*0.000001;
   return (nd); 
 }
  double xznd(double nd,double r0,double Rs)  /******t溫度井筒內(nèi)原油的修正粘度******/
  {
	  double x,A,B,ndp,nd0;
    x=0.15*7.15*(0.0239+0.01638*pow(nd/r0,0.278));
    ndp=nd*pow(10,x);
    A=4.404*pow((Rs+17.7935),-0.515);
    B=3.0355*pow((Rs+26.6905),-0.338);
    nd0=A*pow(ndp,B);
    return (nd0); 
  }             
 double bmzl(double r,double t,double p)  /******t溫度井筒內(nèi)原油的表面張力******/
 { 
	 double xx,yy,zl;
     xx=0.267*(141.5-131.5*r)/r;
	 yy=42.4-0.047*(1.8*t+32)-xx;
	 zl=y*exp(-0.0001015*p);
	 return (zl);
 }
 double wlmz(double Re,double d)   /******水力光滑******/
 {  
	 double aa,ee,xx,yy,zz;
	 double mz;
     aa=0.19;
     ee=(2*a)/d;
	if(Re>3000&&Re<(59.7/pow(ee,8/7))) /******水力光滑******/
      mz=0.3164/pow(Re,0.25);
	if(Re>(59.7/pow(ee,8/7))&&Re<((665-765*log10(ee))/ee))
	{  
      xx=6.8/Re+pow(aa/(3.7*d),1.11);
	  yy=-1.8*log10(xx);
	  mz=1/(yy*yy);
	}
	if(Re>((665-765*log10(ee))/ee))
    { 
      zz=2*log10(3.7*d/aa);
	  mz=1/(zz*zz);
	}
	return (mz);
}
 main()
{ 
  double tf[40],t[40],k[40],tw[40],tr[40],tm[40],xq[40],rq[40],
        ndq[40],bq[40],Ra[40],xh[40];
  double  p0,d1,d2,d3,d4,d5,W,T,fs,cp,kl;/******溫度的參數(shù)*****/
  double Rp,Rs,G1,G2,G,Ap,pl,nd1,nd2;
  double r0[40],r01[40],p[40],pm[40],pj[40],Rs1[40],Rs2[40],Rs0[40],
         B0[40],B01[40],ndz[40],q0[40],V[40],v[40],Re[40],f[40],rm[40];/******單相流的參數(shù)*****/
  double vsg,w,N,N1,E1,L1,L2,L3,L4,a1,b,c,d,e,f0,g0,H,H1,Hfl,Hjx,H1fl,cfl,H1jx,cjx,
         cc,A,B,r12,y0,S1,S,zl1,zl,Rel,Rel1,Rel2,fz,fm;  /******多相流的參數(shù)*****/
  int i;
  FILE *fp1;
  p0=16; d1=0.062; d2=0.073;d3=0.1217;
   d4=0.1397; d5=0.2472; cp=1007; nd2=0.00122;
  fp1=fopen("溫度和壓力.txt","w");  /*****打開文件*****/
  printf("井筒內(nèi)的溫度和壓力分布：\n");
  fprintf(fp1,"井筒內(nèi)的溫度和壓力分布：\n");
  printf("i    h      t          p\n");
  fprintf(fp1,"i    h      t           p\n");
  T=30.87+0.508*C+5;  /*****原油的凝固點*****/
  for(i=1;i<=40;i++)
  {
	W=1000*y*cp1/24+y*GOR*cp2/24;
    W=W*1000/3600;  /*****混合物的水當量*****/
    tf[i]=t0-m*i*50; /*****50*i處的地層溫度*****/
    fs=0.9812*log(1+1.81*2*sqrt(a*s*24*3600)/d5);  /******無因次溫度*****/
  do{  
    k[i]=2.5;
    t[i]=W*m*(1-exp(-k[i]*50*i/W))/k[i]+tf[i];  /*****油氣混合物的溫度*****/
    if(t[i]<T)
     t[i]=(W*m+qi)*(1-exp(-k[i]*50*i/W))/k[i]+tf[i]; /*****電加熱時油氣混合物的溫度*****/
     tw[i]=(xd*tf[i]+t[i]*k[i]*fs*d2/2)/(xd+k[i]*fs*d2/2); /******水泥環(huán)外緣的溫度*****/
     tr[i]=tw[i]+k[i]*(t[i]-tw[i])*log(d5/d4)/(2*PI*xs); /******水泥環(huán)與套管處的溫度*****/
     tm[i]=(tr[i]+tw[i])/2;  /******環(huán)空中平均溫度*****/
     xq[i]=(2.83+0.0073*(tm[i]-50))*0.01; /******環(huán)空中氣體導熱系數(shù)*****/
     rq[i]=1.093-0.0031*(tm[i]-50);  /******環(huán)空中氣體密度*****/ 
     ndq[i]=(19.6+0.05*(tm[i]-50))*1e-6; /******環(huán)空中氣體粘度*****/
     bq[i]=1/(273+tm[i]); /******環(huán)空中氣體膨脹系數(shù)*****/
     Ra[i]=g*bq[i]*rq[i]*rq[i]*cp*(t[i]-tr[i])*pow((d3-d2)/2,3)/(ndq[i]*xq[i]);/******判斷環(huán)空當量導熱系數(shù)的瑞利準數(shù)*****/
     if(Ra[i]<6000)
     xh[i]=xq[i];
     if(Ra[i]>6000&&Ra[i]<2e5)
     xh[i]=0.13*xq[i]*pow(Ra[i],0.25);
     if(Ra[i]>2e5&&Ra[i]<1.1e7)
     xh[i]=0.048*xq[i]*pow(Ra[i],0.333);
     kl=1/(log(d3/d2)/(2*PI*xh[i])+log(d5/d4)/(2*PI*xs)); /******根據(jù)熱阻求得的傳熱系數(shù)*****/
     if(fabs(k[i]-kl)>0.001)
      k[i]=kl;
     break; 
  }while(1);

        /********壓力的計算*******/
      Rp=GOR*r1;  /******產(chǎn)出氣油比*****/
      r0[i]=r1-0.000699*(tm[i]-20);  /******平均相對密度*****/
      r01[i]=1000*r0[i];   /******平均密度*****/
      G1=y*1000/24;        /******日產(chǎn)油量*****/
      G2=y*GOR*r2*1.293/24;   /******日產(chǎn)氣量*****/
      G=G1+G2;       /******日產(chǎn)油氣量*****/
      Ap=0.7854*d1*d1;   /******油管截面積*****/
      w=G/Ap/r01[i]/3600;        /******混合物流速*****/
      nd1=niandu(tm[i]);
      Rs1[i]=pow(10,(1.76875/r0[i]-0.00168*tm[i]));
      Rs2[i]=7.15*Rs1[i]/5.96997;
      Rs0[i]=0.1745*r2*pow(Rs2[i],1/0.83);   /******溶解氣油比*****/
      p[i]=15;
      pm[i]=(p[i]+p0)/2;  /******流段平均壓力*****/
      B01[i]=7.1174+Rs0[i]*pow(r2/r0[i],0.5)+0.4005*tm[i];
      B0[i]=0.972+1.1175*0.001*pow(B01[i],1.175);  /******體積系數(shù)*****/
   if(Rp<=Rs0[i])    /******單相流*****/
   do{ 
	   p[0]=16;
	   V[i]=G1*B0[i]/r01[i];    /******日產(chǎn)油氣總體積*****/
       rm[i]=G/V[i];    /******油氣平均流速*****/
       ndz[i]=xznd(nd1,r0[i],Rs0[i]);   /******調(diào)用函數(shù)求粘度*****/
       q0[i]=G1*B0[i]/(r01[i]*3600);   /******原油體積流量*****/
       v[i]=q0[i]/Ap;    /******流體表觀流速*****/
       Re[i]=r01[i]*d1*v[i]/ndz[i];   /******雷諾數(shù)*****/
       if(Re[i]<3000)
       f[i]=64/Re[i];
	   if(Re[i]>3000)
       f[i]=wlmz(Re[i],d1);
       pj[i]=(f[i]*50*tm[i]*v[i]*v[i]/(2*d1)+rm[i]*50*g)*0.000001;  /******總壓降*****/
       pl=p[i-1]-pj[i];
       if(fabs(p[i]-pl)>0.01);
           p[i]=pl;
          break;  
   }while(1);
     if(Rp>Rs0[i])   /******多相流*****/
      do{ p[0]=16;
          vsg=y*GOR/24/3600;    /******氣相表觀流速*****/ 
          N=w*w/(g*d1);
          E1=G1/(G1+G2);
          L1=316*pow(E1,0.302);
          L2=92.52*0.00001*pow(E1,-2.4684);
          L3=0.1*pow(E1,-1.4516);
          L4=0.5*pow(E1,-6.733);
		  N1=bmzl(r1, tm[i],p[i]);
      if((E1<0.01&&N<L1)||(E1>=0.01&&N<L2))  /******分離流******/
       { 
		 a1=0.98;b=0.4846;c=0.0868;d=0.011;
         e=-3.768;f0=3.539;g0=-1.614;
         H1=a1*pow(E1,b)/pow(N,c);
        if(H1<E1)
		   H1=E1;
           cc=(1-E1)*log(d*pow(E1,e)*pow(N1,f0)*pow(N,g0));
           H=H1*(1+0.33*cc);  }
      if(((E1>=0.01&&E1<0.4)&&(N>L3&&N<L1))||(E1>=0.4&&(N>L3&&N<L4)))  /******間歇流******/
       {  
		  a1=0.845;b=0.5351;c=0.0173;d=2.96;
          e=0.305;f0=-0.4473;g0=0.0978;
          H1=a1*pow(E1,b)/pow(N,c);
           if(H1<E1) H1=E1;
           cc=(1-E1)*log(d*pow(E1,e)*pow(N1,f0)*pow(N,g0));
           H=H1*(1+0.33*cc);  
	  }
       if((E1<0.4&&N>=L1)||(E1>=0.4&&N>L4))    /******分散流******/
      {   
		   a1=1.065;b=0.5929;c=0.0609;
          H=a1*pow(E1,b)/pow(N,c);
	   }
       if(E1>=0.01&&(N>L2&&N<=L3))     /******過渡流******/
      {  
		   H1fl=0.98*pow(E1,0.4846)*pow(N,0.0868);
          if(H1fl<E1) H1fl=E1;
          cfl=(1-E1)*log(0.011*pow(E1,-3.768)*pow(N1,3.539)*pow(N,-1.614));
          Hfl=H1fl*(1+0.33*cfl);
          H1jx=0.845*pow(E1,0.5351)*pow(N,0.0173);
          if(H1jx<E1) H1jx=E1;
          cjx=(1-E1)*log(2.96*pow(E1,0.305)*pow(N1,-0.4473)*pow(N,0.0978));
          Hjx=H1jx*(1+0.33*cjx);
           A=(L3-N)/(L3-L2);
           B=1-A;
           H=A*Hfl+B*Hjx;  
	   }
           r12=(r1*H+r2*1.293*pm[i]*t[i]*(1-H)/(0.899*tm[i]*p[i]))*1000; /******混合物實際密度******/
           y0=E1/(H*H);
           S1=-0.0523+3.18*log(y0)-0.8725*pow(log(y0),2)+0.01853*pow(log(y0),4);
           Rel1=d1*w*(r1*E1+r2*(1-E1));
           Rel2=nd1*E1+nd2*(1-E1);
           Rel=Rel1/Rel2;
           zl1=0.0056+0.5/pow(Rel,0.32);
           zl=zl1*exp(S1);   /******兩相流阻力系數(shù)******/
           fz=r12*g+zl*(G*1000/24/3600)*w/(2*d1*Ap);
           fm=1-r12*w*vsg/pm[i]/100000;
           pj[i]=(fz/fm)*50*0.000001;     /******總壓降*****/
           pl=p[i-1]-pj[i];
        if(fabs(p[i]-pl)>0.01);
           p[i]=pl;
          break;  
	  }while(1);
       printf("%d   %d    %4.2f    %10.4f\n",i,50*i,t[i],p[i]);
       fprintf(fp1,"%d   %d    %4.2f    %10.4f\n",i,50*i,t[i],p[i]);
        }
      fclose(fp1);
       getch();     }