****** TO LOAD CASE 108:TYPE L(N108) *****
    GROUP 1. Run title and other preliminaries
TEXT(2D LAMINAR FLOW OVER A FENCE: N108
TITLE
mesg(PC486/50 time last reported as appx. 16 min
  DISPLAY
  The case considered is 2d laminar incompressible flow over a
  thin fence of height H located in a planar channel with a
  blockage ratio S/H=0.75. The flow Reynolds number is 82.5 based
  on fence height S and inlet bulk velocity Uin. This situation
  has been studied experimentally and numerically by Carvalho et
  al (Appl.Math.Modelling, Vol.11, p23, (1987) ). The boundary
  conditions correspond to an inlet flow of fully-developed
  parabolic velocity profile located at 8 fence heights upstream
  of the fence, and an outlet condition of fixed pressure at 15
  fence heights downstream of the fence, and no-slip conditions
  at the walls. In practice the fence thickness is t/H=0.133,
  but in the calculations it is taken as zero.
 
  The present calculation is performed in all six planes as a
  consistency test. The default calculation is made with the
  linear-upwind scheme (LUS) for momentum, and the option exists
  to make the calculation with the HYBRID scheme. A mesh of 24
  vertical by 48 horizontal cells has been used in the
  calculations.
  ENDDIS
 
  Experiments indicate that a primary recirculation zone develops
  behind the fence with a reattachment length of X1/S=4.4, as
  measured from the fence. An additional separation cell forms on
  the upper wall of the channel with the separation point located
  at X2/S=4.0 and the reattachment point at X3/S=6.0, yielding a
  separation length of DX/S=2.0. The main results may be summarised
  as follows:
                 Data   Hybrid  Linear-Upwind
      X1/S       4.4     5.56       4.96
      X2/S       4.0     4.21       3.85
      X3/S       6.0     6.72       8.62
      DX/S       2.0     2.51       4.77
 
  Although no grid-refinement studies have been performed, the
  results are in reasonable agreement with the data.
  PHOTON USE
  P
 
 
   0.20443E+04 0.15633E+04 CR
  GR OU Z 1;GR Z 1 COL 3;GR X 13 13 Y 1 12 COL 15
  STREAM 2D Z 1 X 4 35 Y 1 M
  -.199E-2 .642E-3 10
  set order
  1
  4
  5
  2
  6
  7
  3
  MAG GR 2
   0.16213E+04 0.17326E+04 CR
  msg Streamlines
  msg press 
  pause
  stream off;grid off 2;redr
  msg vectors
  vec z 1 sh
  msg press  and then  to end
  pause
  ENDUSE
  AUTOPLOT USE
  file
  phi 5
 
  d 1 u1 y 1;d 1 u1 y m;div x 7.5e-3 1 2;shift x -8 1 2
  col3 1;colf 2;level y 0;scale x 0 10;scale y -.3 .25;
  msg horizontal velocity along bottom & top walls behind fence
  msg separation & reattachment points when U1 passes through zero
  msg blue line - U1 along bottom wall
  msg red  line - U1 along top wall
  msg press  to continue
  msg press  to end
  ENDUSE
 
REAL(YG,YFP,YFM,CWID2,DTF,CWID,CDOWN,CUP,CLEN,SFENCE,REYNO)
REAL(UBIN,UIN,SBLOCK)
INTEGER(NVFENC,NVGAP,NHUP,NHDOWN)
CHAR(SCHM,DIRV,DIRH,VELV,VELH,PLANE)
     ** Calculation of domain specifications
CWID=0.01;SBLOCK=0.75;SFENCE=SBLOCK*CWID;CUP=8.*SFENCE
CDOWN=15.*SFENCE;CLEN=CUP+CDOWN;REYNO=82.5
NVFENC=10+2;NVGAP=10+2;NHUP=12;NHDOWN=28;nhdown=36
ENUL=1.8E-5;RHO1=1.2;UBIN=REYNO*ENUL/SFENCE
CWID2=CWID*CWID
MESG( Enter required vertical coordinate X,Y or Z
MESG( Default:  Y
READVDU(DIRV,CHAR,Y)
CASE :DIRV: OF
WHEN X,1
+ NREGX=2;VELV=U1
WHEN Y,1
+ NREGY=2;VELV=V1
WHEN Z,1
+ NREGZ=2;VELV=W1
ENDCASE
MESG( Enter required horizontal coordinate X,Y or Z
MESG( Default:  X
READVDU(DIRH,CHAR,X)
CASE :DIRH: OF
WHEN X,1
+ NREGX=2;VELH=U1
WHEN Y,1
+ NREGY=2;VELH=V1
WHEN Z,1
+ NREGZ=2;VELH=W1
ENDCASE
PLANE=:DIRV::DIRH:
TEXT(2D:PLANE: LAMINAR FLOW OVER A FENCE: N108
    GROUP 3. X-direction grid specification
    GROUP 4. Y-direction grid specification
     ** channel length = 0.762 & channel width = 0.1143
IREG:DIRH:=1;GRDPWR(:DIRH:,NHUP,CUP,-1.2)
IREG:DIRH:=2;GRDPWR(:DIRH:,NHDOWN,CDOWN,1.3)
IREG:DIRV:=1;GRDPWR(:DIRV:,NVFENC,SFENCE,-1.2)
IREG:DIRV:=2;GRDPWR(:DIRV:,NVGAP,CWID-SFENCE,-1.3)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,:VELH:,:VELV:);SOLUTN(P1,Y,Y,Y,P,P,P)
    GROUP 8. Terms (in differential equations) & devices
MESG( Enter required convection scheme
MESG( Default: LUS - Linear upwind for momentum Scheme
MESG( The alternative is:
MESG(  HYB  - Hybrid Differencing Scheme for all variables
READVDU(SCHM,CHAR,LUS)
CASE :SCHM: OF
WHEN HYB,3
+ MESG(Hybrid-differencing scheme
+ DIFCUT=0.5;DTF=CLEN/UBIN
+ LSWEEP=800
WHEN LUS,3
+ MESG(Linear upwind scheme for momentum
+ SCHEME(LUS,:VELH:,:VELV:)
+ DTF=0.06*CLEN/(UBIN*N:DIRH:)
+ LSWEEP=1600
ENDCASE
    GROUP 9. Properties of the medium (or media)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:VELH:)=UBIN;FIINIT(:VELV:)=0.1*UBIN
    GROUP 13. Boundary conditions and special sources
IF(:PLANE:.EQ.YX) THEN
+ CONPOR(FENCE,0.0,EAST,-%1,-%1,-#1,-#1,#1,#1)
+ PATCH(OUTLET,EAST,#NREGX,#NREGX,#1,#NREGY,#1,#1,1,1)
+ WALL (TOP,NORTH,#1,#NREGX,#NREGY,#NREGY,#1,#1,1,1)
+ WALL (BOT,SOUTH,#1,#NREGX,#1,#1,#1,#1,1,1)
  ** Inlet velocity profile: u(y)/uav = {6y(cwid-y)}/cwid**2
     umax/uav=1.5 at y=cwid/2.
+ YFM=0.
DO JJ=1,NY
+ YFP=YFRAC(JJ)*YVLAST;YG=0.5*(YFP+YFM)
+ UIN=6.0*UBIN*YG*(CWID-YG)/CWID2;YFM=YFP
+ PATCH(IN:JJ:,WEST,1,1,JJ,JJ,1,NZ,1,LSTEP)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*UIN)
+ COVAL(IN:JJ:,U1,ONLYMS,UIN);COVAL(IN:JJ:,V1,ONLYMS,0.0)
ENDDO
ENDIF
IF(:PLANE:.EQ.XY) THEN
+ CONPOR(FENCE,0.0,NORTH,-#1,-#1,-%1,-%1,#1,#1)
+ PATCH(OUTLET,NORTH,#1,#NREGX,#NREGY,#NREGY,#1,#1,1,1)
+ WALL (TOP,EAST,#NREGX,#NREGX,#1,#NREGY,#1,#1,1,1)
+ WALL (BOT,WEST,#1,#1,#1,#NREGY,#1,#1,1,1)
+ YFM=0.
DO JJ=1,NX
+ YFP=XFRAC(JJ)*XULAST;YG=0.5*(YFP+YFM)
+ UIN=6.0*UBIN*YG*(CWID-YG)/CWID2;YFM=YFP
+ PATCH(IN:JJ:,SOUTH,JJ,JJ,1,1,1,NZ,1,LSTEP)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*UIN)
+ COVAL(IN:JJ:,U1,ONLYMS,0.0);COVAL(IN:JJ:,V1,ONLYMS,UIN)
ENDDO
ENDIF
IF(:PLANE:.EQ.YZ) THEN
+ CONPOR(FENCE,0.0,HIGH,#1,#1,-#1,-#1,-%1,-%1)
+ PATCH(OUTLET,HIGH,1,1,#1,#NREGY,#NREGZ,#NREGZ,1,1)
+ WALL (TOP,NORTH,1,1,#NREGY,#NREGY,#1,#NREGZ,1,1)
+ WALL (BOT,SOUTH,1,1,#1,#1,#1,#NREGZ,1,1)
+ YFM=0.
DO JJ=1,NY
+ YFP=YFRAC(JJ)*YVLAST;YG=0.5*(YFP+YFM)
+ UIN=6.0*UBIN*YG*(CWID-YG)/CWID2;YFM=YFP
+ PATCH(IN:JJ:,LOW,1,NX,JJ,JJ,1,1,1,LSTEP)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*UIN)
+ COVAL(IN:JJ:,W1,ONLYMS,UIN);COVAL(IN:JJ:,V1,ONLYMS,0.0)
ENDDO
ENDIF
IF(:PLANE:.EQ.ZY) THEN
+ CONPOR(FENCE,0.0,NORTH,#1,#1,-%1,-%1,-#1,-#1)
+ PATCH(OUTLET,NORTH,1,1,#NREGY,#NREGY,#1,#NREGZ,1,1)
+ WALL (TOP,HIGH,1,1,#1,#NREGY,#NREGZ,#NREGZ,1,1)
+ WALL (BOT,LOW,1,1,#1,#NREGY,#1,#1,1,1)
+ YFM=0.
DO JJ=1,NZ
+ YFP=ZFRAC(JJ)*ZWLAST;YG=0.5*(YFP+YFM)
+ UIN=6.0*UBIN*YG*(CWID-YG)/CWID2;YFM=YFP
+ PATCH(IN:JJ:,SOUTH,1,1,1,1,JJ,JJ,1,LSTEP)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*UIN)
+ COVAL(IN:JJ:,V1,ONLYMS,UIN);COVAL(IN:JJ:,W1,ONLYMS,0.0)
ENDDO
ENDIF
IF(:PLANE:.EQ.XZ) THEN
+ CONPOR(FENCE,0.0,HIGH,-#1,-#1,#1,#1,-%1,-%1)
+ PATCH(OUTLET,HIGH,#1,#NREGX,1,1,#NREGZ,#NREGZ,1,1)
+ WALL (TOP,EAST,#NREGX,#NREGX,1,1,#1,#NREGZ,1,1)
+ WALL (BOT,WEST,#1,#1,1,1,#1,#NREGZ,1,1)
+ YFM=0.
DO JJ=1,NX
+ YFP=XFRAC(JJ)*XULAST;YG=0.5*(YFP+YFM)
+ UIN=6.0*UBIN*YG*(CWID-YG)/CWID2;YFM=YFP
+ PATCH(IN:JJ:,LOW,JJ,JJ,1,NY,1,1,1,LSTEP)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*UIN)
+ COVAL(IN:JJ:,W1,ONLYMS,UIN);COVAL(IN:JJ:,U1,ONLYMS,0.0)
ENDDO
ENDIF
IF(:PLANE:.EQ.ZX) THEN
+ CONPOR(FENCE,0.0,EAST,-%1,-%1,#1,#1,-#1,-#1)
+ PATCH(OUTLET,EAST,#NREGX,#NREGX,1,1,#1,#NREGZ,1,1)
+ WALL (TOP,HIGH,#1,#NREGX,1,1,#NREGZ,#NREGZ,1,1)
+ WALL (BOT,LOW,#1,#NREGX,1,1,#1,#1,1,1)
+ YFM=0.
DO JJ=1,NZ
+ YFP=ZFRAC(JJ)*ZWLAST;YG=0.5*(YFP+YFM)
+ UIN=6.0*UBIN*YG*(CWID-YG)/CWID2;YFM=YFP
+ PATCH(IN:JJ:,WEST,1,1,1,NY,JJ,JJ,1,LSTEP)
+ COVAL(IN:JJ:,P1,FIXFLU,RHO1*UIN)
+ COVAL(IN:JJ:,U1,ONLYMS,UIN);COVAL(IN:JJ:,W1,ONLYMS,0.0)
ENDDO
ENDIF
COVAL(OUTLET,P1,1.E3,0.0)
COVAL(OUTLET,:VELH:,ONLYMS,0.0);COVAL(OUTLET,:VELV:,ONLYMS,0.0)
    GROUP 17. Under-relaxation devices
RELAX(:VELH:,FALSDT,DTF);RELAX(:VELV:,FALSDT,DTF)
I:DIRV:MON=NVFENC-2;I:DIRH:MON=NHUP+12;NPRMON=100
    GROUP 23. Field print-out and plot control
TSTSWP=-1;ITABL=3;NPLT=20;IPLTL=LSWEEP
  NSAVE=PH:PLANE: