TALK=T;RUN( 1, 1)
  
  
  

 
TEXT(S-A 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 high
    Reynolds form of the Spalart-Allmaras turbulence model. A
    similar case using the standard k-e turbulence model is 
    available in Library Case b583.
	
    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.
	
    Note that this case must be run in the command line and
    not through the VR Viewer.
  ENDDIS
 
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
DIAM=0.889
 
L1DD=4.0
L3DD=4.0
 
LENG1= l1dd*DIAM
LENG2= -RCURV*ANGLE
LENG3= l3dd*DIAM
 
DIAMEB=0.874
 
 
   **Grid
POWZ1=0.8; POWZ2=1.0; POWZ3=1.1
INTEGER(LASTK,FIRSTK,NZ1,NZ2,NZ3,NY1,NY2)
REAL(POWY1,POWY2)
NY1=16;NY2=2;POWY1=0.7;POWY2=1.0
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
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:)
 
    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:)
 
CONPOR(BLOK1,0.0,CELL,-1,-NX,-#2,-#2,-#2,-#2)
 
    **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)
 
 
qflow=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(SPALART-ALLMARAS)
STORE(ENUT,YPLS,UC1,VC1,WC1,PRPS)
FIINIT(ENTI)= 0.012474
FIINIT(U1)=0.0; FIINIT(V1)=0.0; FIINIT(W1)=0.0
    ** Inlet
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,ENTI,ONLYMS,0.012474)
    ** Outlet
PATCH(OUT,HIGH,1,NX,1,NY,NZ,NZ,1,1);COVAL(OUT,P1,1.0E3,0)
    ** Pipe wall
PATCH(WN1,NWALL,1,NX,NY,NY,#1,#1,1,1)
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,U1,GRND2,0.0);COVAL(WN3,W1,GRND2,0.0)
 
    ** Relaxation parameters
DTHYD=5.*LENG1/(NZ1*WIN)
  RELAX(U1,FALSDT,DTHYD);RELAX(V1,FALSDT,DTHYD)
  RELAX(W1,FALSDT,DTHYD)
RELAX(ENTI,LINRLX,0.5)
    ** 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=500
   ** 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
dref=diameb
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
   ** 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,ENTI ,C,MINMOD)
EGWF=F
WALPRN=T
 
DISTIL=T
EX(P1  )=   3.656E+02
EX(U1  )=   3.369E-01
EX(V1  )=   1.703E-01
EX(W1  )=   1.319E+01
EX(UC1 )=   1.605E-01
EX(VC1 )=   7.347E+00
EX(WC1 )=   7.499E+00
EX(PRPS)=   9.573E-01
EX(YPLS)=   4.751E+02
EX(ENUT)=   4.795E-03
EX(DWDY)=   9.573E-11
EX(DWDX)=   9.573E-11
EX(DVDZ)=   9.573E-11
EX(DVDX)=   9.573E-11
EX(DUDZ)=   9.573E-11
EX(DUDY)=   9.573E-11
EX(LTLS)=   2.516E-02
EX(WDIS)=   1.309E-01
EX(ENTI)=   4.795E-03
EX(WCRT)=   7.499E+00
EX(VCRT)=   7.347E+00
EX(UCRT)=   1.605E-01
EX(VPOR)=   9.573E-01
 LIBREF = 0
STOP