TEXT(3D TURBULENT FLOW IN AN ELBOW METER
TITLE
  DISPLAY
    The case considered is the 3D steady, incompressible,
    isothermal turbulent flow of natural gas at 44bar in a circular
    elbow meter. The pipe inlet diameter is 0.889m, the pipe bend
    is 0.874m diameter and the mean radius of curvature of the bend
    is 1.3145m. The computational domain is extended 4 inlet
    diameters upstream and downstream of the elbow. The volumetric
    flow rate is 8.57 m3/s, the gas compressibility factor is 0.824,
    and the inlet Reynolds number is 3.52E7. The elbow meter is a
    flow-measurement device which exploits the fact that the flow
    rate through the bend is proportional to radial pressure
    difference between the inside and outside of the bend. The
    calculations are carried out with the GCV solver and the k-e
    turbulence model. The computed radial pressure difference at the
    symmetry plane and 22.5o into the elbow is 54mbar, which is in
    excellent agreement with that given by the gas-industry
    correlation of flow rate versus pressure difference.
  ENDDIS
  PHOTON USE
    p
  phi
 
  xyz
   0.00000E+00 0.00000E+00 CR
  UP Z
  VIEW X
  CON P1 X 1 1 Y 1 M Z 10 32 FI;0.01
  GR OU X 1
  CON P1 X m m Y 1 M Z 10 32 FI;0.01
  GR OU X m
  gr z 23
  pause
  mag gr 3
   0.78183E+03 0.26641E+04 CR
  REDR
  pause
  cl;vec x 1 1 y 1 m z 10 32 sh
  vec x m m y 1 m z 10 32 sh
  redr
 
 
  ENDUSE
 
  ** BLOCK=T introduces sudden contraction in elbow
  ** PLUG =T uses uniform inlet profiles
          =F uses fully-develped inlet velocity profile
BOOLEAN(BLOCK,PLUG);BLOCK=T;PLUG=T
REAL(DIAM,LENG1,LENG2,LENG3,ANGLE,POWY,POWZ1,POWZ2,POWZ3)
REAL(L1DD,L3DD,RCURV,PI,DIAMEB,RADIN,RADEB)
INTEGER(KFF,KLL)
PI=3.14159
   ** Pipe parameters
RCURV= 1.3145
ANGLE=-0.5*PI
mesg(Enter pipe inlet diameter (default=0.889m)
readvdu(diam,real,0.889)
 
mesg(
mesg(Enter upstream pipe length (default=4 diameters)
readvdu(l1dd,real,4.0)
mesg(
mesg(Enter downstream pipe length (default=4 diameters)
readvdu(l3dd,real,4.0)
 
LENG1= l1dd*DIAM
LENG2= -RCURV*ANGLE
LENG3= l3dd*DIAM
 
IF(BLOCK) THEN
+ mesg(Enter pipe elbow diameter (default=0.874m)
+ readvdu(diameb,real,0.874)
ENDIF
 
   **Grid
POWZ1=0.8; POWZ2=1.0; POWZ3=1.1
INTEGER(LASTK,FIRSTK,NZ1,NZ2,NZ3,NY1,NY2)
REAL(POWY1,POWY2)
IF(BLOCK) THEN
+ NY1=16;NY2=2;POWY1=0.7;POWY2=1.0
ELSE
+ NY=15;POWY=1.01
ENDIF
NX =15
NZ1=12; NZ2=15; NZ3=12
    GROUP 3. X-direction grid specification
CARTES=F; XULAST=PI
GRDPWR(X,NX,XULAST,1.0)
    GROUP 4. Y-direction grid specification
IF(BLOCK) THEN
+ RADIN =0.5*DIAM;RADEB=0.5*DIAMEB
+ NREGY=2
+ IREGY=1; GRDPWR(Y,NY1,RADEB,:POWY1:)
+ IREGY=2; GRDPWR(Y,NY2,RADIN-RADEB,:POWY2:)
ELSE
+ YVLAST=DIAM*0.5
+ GRDPWR(Y,NY,YVLAST,:POWY:)
ENDIF
    GROUP 5. Z-direction grid specification
ZWLAST=0.1
GRDPWR(Z,10,ZWLAST,1.0)
NREGZ = 3
IREGZ=1; GRDPWR(Z,NZ1,LENG1,:POWZ1:)
IREGZ=2; GRDPWR(Z,NZ2,LENG2,1.0)
IREGZ=3; GRDPWR(Z,NZ3,LENG3,1.0)
    GROUP 6. Body-fitted coordinates or grid distortion
BFC= T
KFF=NZ1+1; KLL=NZ1+NZ2+1
  **
GSET(C,K:KLL :,F,K:KFF:,1,NX,1,NY,RX,ANGLE,RCURV,LENG1,INC,:POWZ2:)
GSET(C,K:NZ+1:,F,K:KLL:,1,NX,1,NY,+,0.0,LENG3,0.0,INC,:POWZ3:)
 
IF(BLOCK) THEN
+ CONPOR(BLOK1,0.0,CELL,-1,-NX,-#2,-#2,-#2,-#2)
ENDIF
    **Physical characteristics
real(qflow,ain,kein,epin,win,rey,dthyd,zfact,gascon,pin,tin)
real(wmax,an,ran,wcur,xcur,ycur,zcur,rcur,xcen,ycen,zcen)
pin=44.e5;tin=(18.7+273);gascon=8314.3/16.;zfact=0.8235
integer(jj1)
enul= 3.47E-7;rho1=pin/(zfact*gascon*tin)
 
mesg(Enter volumetric flow rate (default=8.57 m^3/s)
readvdu(qflow,real,8.57)
 
ain = 0.25*pi*diam*diam ; win  = qflow/ain
rey = win*diam/enul
 
real(fric,aa,bb)
  ** use Karman-Nikuradse correlation to estimate f
fric=0.25/(1.82*log10(rey)-1.64)**2
do jj=1,20
+ bb=sqrt(fric);bb=rey*bb;aa=log10(bb)
+ bb=4.0*aa-0.4;fric=1./(bb*bb)
enddo
  ** use Hinze data to estimate average turbulence levels
kein = fric*win*win
epin = 0.1643*kein**1.5/(0.045*diam)
 
    **Flow settings.
SOLVE(P1,U1,V1,W1);TURMOD(KEMODL)
FIINIT(KE)=KEIN; FIINIT(EP)=EPIN
FIINIT(U1)=0.0; FIINIT(V1)=0.0; FIINIT(W1)=0.0
    ** Inlet
IF(PLUG) THEN
+ PATCH(IN,LOW,1,NX,1,NY,1,1,1,1)
+ COVAL(IN,P1,FIXFLU,RHO1*WIN);COVAL(IN,W1,ONLYMS,WIN)
+ COVAL(IN,U1,ONLYMS,0.0);COVAL(IN,V1,ONLYMS,0.0 )
+ COVAL(IN,KE,ONLYMS,KEIN);COVAL(IN,EP,ONLYMS,EPIN)
ELSE
+ AN=1./(4.*FRIC)**0.5
+ RAN=1./AN
+ WMAX=WIN*0.5*(AN+1.)*(2*AN+1.)*RAN*RAN
WMAX
+ XCEN= 0.5*(XC(1,1,1)+XC(2,1,1))
+ YCEN= 0.5*(YC(1,1,1)+YC(2,1,1))
+ ZCEN= 0.5*(ZC(1,1,1)+ZC(2,1,1))
DO JJ=1,NY
+ PATCH(IN:JJ:,LOW,1,NX,JJ,JJ,1,1,1,1)
+ JJ1 = JJ+1
+ XCUR= 0.25*(XC(1,JJ,1)+XC(2,JJ,1)+XC(1,JJ1,1)+XC(2,JJ1,1))
+ YCUR= 0.25*(YC(1,JJ,1)+YC(2,JJ,1)+YC(1,JJ1,1)+YC(2,JJ1,1))
+ ZCUR= 0.25*(ZC(1,JJ,1)+ZC(2,JJ,1)+ZC(1,JJ1,1)+ZC(2,JJ1,1))
+ XCUR= XCUR-XCEN ; YCUR= YCUR-YCEN ; ZCUR= ZCUR-ZCEN
+ RCUR= XCUR*XCUR+YCUR*YCUR+ZCUR*ZCUR; RCUR=SQRT(RCUR)
+ WCUR= WMAX*(1.0-RCUR/RADIN)**RAN
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*WCUR)
+ COVAL(IN:JJ:,W1,ONLYMS,WCUR)
+ COVAL(IN:JJ:,U1,ONLYMS,0.0);COVAL(IN:JJ:,V1,ONLYMS,0.0)
+ COVAL(IN:JJ:,KE,ONLYMS,KEIN);COVAL(IN:JJ:,EP,ONLYMS,EPIN)
ENDDO
ENDIF
    ** Outlet
PATCH(OUT,HIGH,1,NX,1,NY,NZ,NZ,1,1);COVAL(OUT,P1,1.0E3,0)
    ** Pipe wall
IF(BLOCK) THEN
+ PATCH(WN1,NWALL,1,NX,NY,NY,#1,#1,1,1)
+ COVAL(WN1,KE,GRND2,GRND2);COVAL(WN1,EP,GRND2,GRND2)
+ COVAL(WN1,U1,GRND2,0.0);COVAL(WN1,W1,GRND2,0.0)
 
+ PATCH(WN3,NWALL,1,NX,NY,NY,#3,#3,1,1)
+ COVAL(WN3,KE,GRND2,GRND2);COVAL(WN3,EP,GRND2,GRND2)
+ COVAL(WN3,U1,GRND2,0.0);COVAL(WN3,W1,GRND2,0.0)
ELSE
+ PATCH(WN1,NWALL,1,NX,NY,NY,1,NZ,1,1)
+ COVAL(WN1,KE,GRND2,GRND2);COVAL(WN1,EP,GRND2,GRND2)
+ COVAL(WN1,U1,GRND2,0.0);COVAL(WN1,W1,GRND2,0.0)
ENDIF
    ** Relaxation parameters
DTHYD=5.*LENG1/(NZ1*WIN)
  RELAX(U1,FALSDT,DTHYD);RELAX(V1,FALSDT,DTHYD)
  RELAX(W1,FALSDT,DTHYD)
RELAX(KE,LINRLX,0.5);RELAX(EP,LINRLX,0.5);KELIN = 3
    ** Output.
OUTPUT(P1,Y,N,P,Y,Y,Y)
OUTPUT(U1,Y,N,N,Y,Y,Y);OUTPUT(V1,Y,N,N,Y,Y,Y);OUTPUT(W1,Y,N,N,Y,Y,Y)
IXMON =NX/2; IYMON=NY/2; IZMON=NZ/2
ITABL=3;NPLT=10;TSTSWP= -1
LSWEEP=300
   ** pressure difference expected from gas-industry
      correlation: q=0.948*sqrt(0.5*Rc/De)*ae*sqrt(2.*dp/rho)
      where ae=pi*De**2/4.
real(dp,ae,cd,dref);cd=0.948
if(block) then
+ dref=diameb
else
+ dref=diam
endif
ae=0.25*pi*dref*dref
dp=(qflow/(ae*cd))**2;dp=dp*dref*rho1/rcurv
MESG(Expected pressure difference across elbow
dp
 
   ** use arithmetic averaging for momentum diffusion
SOLUTN(U1,P,P,P,P,P,N);SOLUTN(V1,P,P,P,P,P,N)
SOLUTN(W1,P,P,P,P,P,N)
   ** activate gcv solver
   
GCV=T
IF(GCV) THEN
   ** use minmod discretisation scheme for convection
+ SPEDAT(SET,GCVSCH,UC1,C,MINMOD)
+ SPEDAT(SET,GCVSCH,VC1,C,MINMOD)
+ SPEDAT(SET,GCVSCH,WC1,C,MINMOD)
+ SPEDAT(SET,GCVSCH,KE ,C,MINMOD)
+ SPEDAT(SET,GCVSCH,EP ,C,MINMOD)
ENDIF
EGWF=F
WALPRN=T;DISTIL=T
IZPRF=23;IZPRL=23;NZPRIN=1;NYPRIN=1;NXPRIN=1
STORE(UC1,VC1,WC1,PRPS)
EX(P1  )=4.092E+02;EX(U1  )=3.443E-01;EX(V1  )=1.972E-01
EX(W1  )=1.327E+01;EX(KE  )=5.664E-01;EX(EP  )=8.216E+00
EX(UC1 )=1.459E-01;EX(VC1 )=7.389E+00;EX(WC1 )=7.558E+00
EX(EPKE)=1.215E+01;EX(WCRT)=7.558E+00;EX(VCRT)=7.389E+00
EX(UCRT)=1.459E-01;EX(VPOR)=9.573E-01;EX(PRPS)=9.573E-01