home *** CD-ROM | disk | FTP | other *** search
/ Hall of Fame / HallofFameCDROM.cdr / oilfield / gasper.lzh / GASPER.FOR < prev    next >
Text File  |  1980-01-02  |  8KB  |  253 lines
  1. c*****************************************************************
  2. c          This program simulates gas well performance
  3. c    It computes either static of flowing bottom
  4. c    hole pressure for both sweet and/or sour gas
  5. c    wells 
  6.     character yes*5,no*5,ans*5,what*5,fname*15
  7.     real h(6),w(6),mu,pt
  8.     common theta,tem,w,pt,co2,h2s
  9.      &  ,h,q,tt,beta,e,dia,mu,sg,dept,ffact,zf
  10. c    read data for Gaussian six points quadrature
  11.     write(*,*)'PLEASE GIVE YOUR NAME IN QUOTES'
  12.         write(*,*)'LIKE '' YOUR NAME'', AND HIT ''ENTER'''
  13.     read(*,24)fname
  14.     write(*,*)' HENRY WELCOMES YOU,  ',fname
  15.     write(*,*)'HAVE YOU LOOK AT THE README FILE'
  16.     write(*,*)'********************************'
  17.     write(*,*)' JUST ANSWER YES OR NO THEN HIT ''ENTER'''
  18.     write(*,*)'********************************'
  19.     read(*,24)what
  20.     if(what.eq.'no')then
  21.     write(*,*)' PLEASE DO BEFORE RUNNING THIS PROGRAM'
  22.     write(*,*)'********************************'
  23.     write(*,*)'TO READ, AT THE PROMPT,ENTER '' TYPE README'''
  24.     stop
  25.     endif
  26.     write(*,1)
  27. 1    format(4x,'THIS PROGRAM COMPUTES STATIC OR FLOWING',//,
  28.      &  4x,'BOTTOM HOLE PRESSURES IN SWEET AND/OR SOUR',//,
  29.      &  4x,'GAS WELLS, TO USE, JUST ANSWER THE QUESTIONS',//
  30.      &  ,4x,'PLEASE DIRECT ALL INQUIRIES TO',//,4x,
  31.      &  ' HENRY A. OHEN , SCHOOL OF PETROLEUM',
  32.      &  //,4x,'& GEOLOGICAL ENGINEERING, UNIVERSITY',//,4x
  33.      &  , 'OF OKLAHOMA, NORMAN OK 73019',//)
  34.     open(10,file='gaus.dat')
  35.     do 2 i=1,6
  36.     read(10,*)h(i),w(i)
  37.  2    continue
  38.     imax=100
  39.     write(*,3)
  40.  3    format(4x,'IS YOUR GAS SWEET',//,4x,
  41.      &  ' ANSWER YES OR NO AND HIT ''ENTER''',//) 
  42.     read(*,24)ans
  43.     if(ans.eq.'no')then
  44.     write(*,4)
  45.  4    format(4x,' YOU MUST KNOW THE MOLE PERCENT OF CARBON DIOXIDE',
  46.      &  //,4x,' AND THE MOLE PERCENT OF HYDROGEN SULPHIDE',//,
  47.      & 'DO YOU HAVE THIS INFORMATION? YES OR NO AND HIT ''ENTER''',//)
  48.     read(*,24)ans
  49.     if(ans.eq.'no')then
  50.     write(*,5)
  51.  5    format(4x,' YOU ARE NOT READY,GET THE REQUIRED INFORMATION'
  52.      &  ,//)
  53.     stop
  54.     else
  55.     write(*,6)
  56. 24      format(a9)
  57.  6    format(4x,'GIVE MOLE FRACTION OF CO2 AND H2S IN THAT ORDER',/
  58.      & ,4x,' HIT ''ENTER'' FOR EACH ENRTY')
  59.     read(*,*)co2,h2s
  60.     endif
  61.     else
  62.     co2=0
  63.     h2s=0
  64.     endif
  65.     write(*,7)
  66.     read(*,*)dept
  67.     write(*,9)
  68.     read(*,*)tt
  69.     write(*,10)
  70.     read(*,*)pt
  71.     write(*,12)
  72.     read(*,*)dia
  73.     write(*,14)
  74.     read(*,*)beta
  75.     write(*,16)
  76.     read(*,*)theta
  77.     write(*,17)
  78.     read(*,*)sg
  79.     write(*,19)
  80.     read(*,*)tol
  81.     write(*,20)
  82.     read(*,*)q
  83.  7    format(4x,'GIVE WELL DEPTH(ft), HIT ''ENTER''')
  84.  9    format(4x,'TUBING HEAD TEMP.(of), HIT ''ENTER''')
  85.  10     format(4x,'GIVE PRESURE (psi), HIT ''ENTER''')
  86.  12     format(4x,'PIPE DIAMETER (inches), HIT ENTER''')
  87.  14     format(4x,'GIVE TEMPERATURE GRADIENT (of/ft), HIT ''ENTER''')
  88.  16     format(4x,'GIVE HORIZONTAL INCLINATION (degree), HIT ''ENTER''')
  89.  17     format(4x,'GIVE GAS GRAVITY HIT ''ENTER''')
  90.  19     format(4x,'GIVE TOLERANCE,1.0e-03,IS FINE, HIT ''ENTER''') 
  91.  20     format(4x, 'GIVE GAS FLOW RATE (MMscfd) , HIT ''ENTER''')
  92.  
  93.     theta=theta/180.0
  94.     e=0.0006
  95.     xguess=pt+0.3*pt
  96.  15    call fcn(xguess,f,fp)
  97.     xx=xguess-f/fp
  98.     if(ij.gt.imax)then
  99.     write(*,*)'ITERATION EXCEEDED IN MAIN LOOP, INCREASE TOL'
  100.     stop
  101.     endif
  102.     ch=(xx-xguess)
  103.     if(abs(ch).gt.tol)then
  104.     xguess=xx
  105.     ij=ij+1
  106.     go to 15
  107.     endif
  108.     pw=xx
  109.     open (13,file='res')
  110.     if(q.eq.0)then
  111.     write(13,*)'    ************************'
  112.     write(13,*)'    THIS IS A STATIC WELL'
  113.     write(13,*)'    ************************'
  114.     else
  115.     write(13,*)'    ************************'
  116.     write(13,*)'    THIS IS A FLOWING WELL'
  117.     write(13,*)'    ************************'
  118.     endif
  119.     write(13,31)zf,mu,ffact,pw
  120.  31    format(4x,'Z-FACTOR IS',4x,f7.4,//,
  121.      &  4x,'GAS VISCOSITY AT BOTTOM IS',4x,f7.4,//,
  122.      &   //,4x,'MOODY FRICTION FACTOR IS'
  123.      &  ,4x,f7.4,//4x,
  124.      &  'PRESURE AT GIVEN DEPT IS',2x,f13.5,//)
  125.     write(*,8)
  126.  8    format(//,2x,'  ****************************************')
  127.     write(*,*)'     FIND YOUR RESULTS IN A FILE CALLED ''res'''
  128.     write(*,*)'     ****************************************'
  129.     write(*,*)'     TO READ ENTER ''TYPE RES'''
  130.     write(*,33)fname
  131.     write(*,*)'****************************'
  132. 33    format(//,4x,'HAVE A GOOD DAY,  ',a15,)
  133.     end
  134. c******************************************************
  135.     subroutine fcn(x,f,fp)
  136.     real h(6),w(6),x,mu,integral
  137.     common theta,tem,w,pt,co2,h2s,
  138.      &  h,q,tt,beta,e,dia,mu,sg,dept,ffact,zf
  139.     tem=tt+beta*dept+460
  140.     call solve(x,integral,fp)
  141.     f=1000*dept-integral
  142.     return
  143.     end
  144.     
  145. c***************************************************
  146.     subroutine solve(p,res,fp)
  147.     real t,res,mu,dia,sg,h(6),w(6)
  148.     common theta,tem,w,pt,co2,h2s
  149.      &  ,h,q,tt,beta,e,dia,mu,sg,dept,ffact,zf
  150.     shc=co2*h2s
  151.     rl=dept/(sin(theta))
  152.     ee=120*(shc**0.9-shc**1.6)+15*(h2s**0.5-h2s**4.0)
  153.     tc=170.491+307.344*sg-ee
  154.     pc=709.604-58.715*sg
  155.     pc=pc*tc/((tc+ee)+h2s*(1.0-h2s)*ee)
  156.     ttc=tem/tc
  157.     pr=p/pc
  158.     call zfact(ttc,pr,z)
  159.     zf=z
  160.     call gvis(ttc,tc,pr,z,pc,sg,fm)
  161.     mu=fm
  162.     if(q.ne.0)then
  163.     ren=20000*q*sg/mu/dia
  164.     call fact(ren,fr)
  165.     ffact=fr
  166.     endif
  167.     a1=(53.34*p/sg/tem/zf+111.1*q*q/(dia**4)/p)
  168.     a2=(2.6665*fr*q*q/4)/(dia**5)
  169.     a3=p*p/tem/tem/1000/zf/zf
  170.     fp=-a1/(a2+a3)
  171.     call gaus(fr,z,p,res)
  172.     return
  173.     end
  174.  
  175.     subroutine fact(ren,f)
  176. c******************************************************
  177. c     Subroutine calculates the friction factor by solving
  178. c     The Colebrooks and White Implicit Equation for friction factor
  179. c     using the Newton Raphson's method
  180. c*****************************************************
  181.     real f,ren,b,fn,ff,h(6),w(6)
  182.     common theta,tem,w,pt,co2,h2s
  183.      &  ,h,q,tt,beta,e,dia,mu,sg,dept,ffact,zf
  184. c     Assume Laminar Flow as first approximation
  185.     fo=64/ren
  186.     do 2 i=1,50
  187. c     Start Newton Iteration
  188.     ff=1.74-0.8686*alog(2*e/dia+18.7/sqrt(fo)/ren)-1.0/sqrt(fo)
  189.     aa=(2*e/dia+18.7/ren/sqrt(fo))
  190.     fn=0.5*(1.0+16.24283/aa/ren)/(fo**1.5)
  191.         f=fo-ff/fn
  192.     if(abs((f-fo)/fo).le.1e-05)go to 5
  193.     fo=f
  194.  2    continue
  195.     write(*,*)'iteration failed in ffact after ',i,'times'
  196.     stop
  197.  5    return
  198.     end
  199. c************************************************************
  200.     subroutine gaus(fr,z,xy,res)
  201. c************************************************************
  202.     real y(6),h(6),w(6),length
  203.     common theta,tem,w,pt,co2,h2s
  204.      &  , h,q,tt,beta,e,dia,mu,sg,dept,ffact,zf
  205.     a=pt
  206.     b=xy
  207.     do 20 i=1,6
  208.     y(i)=(h(i)+1.0)*(b-a)/2.0+a
  209.  20    continue
  210.     sum=0
  211.     length=dept/sin(theta)
  212.     do 30 i=1,6
  213.     a1=y(i)*53.34/sg/tem/z
  214.     a2=111.1*q*q/(dia**4)/y(i)
  215.     a3=2.6665*fr*q*q/4/(dia**5)
  216.     a4=y(i)*y(i)/tem/z/z/tem/1000
  217.     g=(a1+a2)/(a3+a4)
  218.     sum=sum+w(i)*g*0.5*(b-a)
  219.  30    continue
  220.     res=sum
  221.     return
  222.     end
  223.  
  224.         subroutine zfact(tr,pr,z)
  225. c       this routine calculates z-factor using  Standing modification
  226. c       to the Brill and Beggs correlation for curve-fitting the Stand-
  227. c       ing and Katz reduced pressure-reduced temperature z-factor chart
  228. c
  229.         a=1.39*(tr-0.92)**0.5-0.36*tr-0.101
  230.         b=(0.62-0.23*tr)*pr
  231.         c=(0.066/(tr-0.86)-0.037)*pr**2
  232.         d=(0.32/(10.**(9.*(tr-1.))))*pr**6
  233.         e=b+c+d
  234.         f=(0.132-0.32*alog10(tr))
  235.         g=10.**(0.3106-0.49*tr+0.1824*tr**2)
  236. c
  237.         z=a+(1.-a)*exp(-e)+f*pr**g
  238.         return
  239.         end
  240. c**********************************************************
  241.  
  242.     subroutine gvis(tr,tc,pr,z,pc,sg,fm)
  243. c     Subroutine computes gas viscosity using the 
  244. c    modified Lee et al correlation
  245.     dg=(2.7*pr*pc*sg)/(z*tr*tc)/(16.10846e03)
  246.     ak=((9.4+0.58*sg)*(tr*pc)**(1.5))/(209.0+
  247.      &  551*sg+tr*tc)
  248.     xx=3.5+986/tc/tr+0.29*sg
  249.     yy=2.4-0.2*xx
  250.     fm=ak*exp(xx*dg**yy)*1e-04
  251.     return
  252.     end
  253.