PHOTON USE
  p
 
 
 
 
  gr x 1
  msg Grid
  vec x 1 sh
  MSG Velocity vectors
  msg Press return to redraw
  pause
  gr cl; gr ou x 1; red
  msg Press return to plot pressure contours
  pause
  cont p1 x 1 fil;.01
  msg Type e to End
  ENDUSE
 
    GROUP 1. Run title
TEXT(PLANE CHANNEL + SMOOTH EXPANSION.:  B533
TITLE
  DISPLAY
    The problem considered is that defined by the International
  Association for Hydraulic Research Working Group on Refined
  Modelling of Flows for their 5th meeting. The report of the
  workshop is presented as a paper entitled 'Laminar flow in a
  complex geometry: a comparison', Int. j. numer. methods fluids,
  vol.5, 667-683(1985).
    The main features of the problem are:
  * laminar flow, Re=10.0 ( and 100.0);
  * the south boundary of the channel is given by the following
     analytical expression :
           y(z) = [tanh(2-30*z/Re)-tanh(2)]/2
           for  0 < z< Re/3
  * the north boundary is a symmetry plane located at y=1;
  * a parabolic profile for the axial velocity is specified at the
    inlet,
           w=3*(y-y*y/2) for z=0 and 0 < y< 1,
           v=0.
  ENDDIS
  Special variables introduced for this problem are:-
    RE........Value of the Reynolds number
    PRESSO....Outlet uniform relative pressure
    CT........Tanh(2)
    POWERZ....Power law in the Z direction
    POWERY....Power law in the Y direction
REAL(RE,CT,A,B,VV,UU,ZZ,POWERZ,POWERY,C,DD,Y0,Y1,YY)
RE=10.0;CT=0.9640275
 
    GROUP 4. Y-direction grid specification
NY=21;POWERY=1.3
    GROUP 5. Z-direction grid specification
NZ=21;POWERZ=1.5;ZWLAST=RE/3.
    GROUP 6. Body-fitted coordinates or grid distortion
BFC=T; NONORT=T; SYMBFC=T
   **The next three YCs are kept at 0.0 to ensure that the w
     resolutes for iz=1 are parallel to the axis. This simplifies
     the specification of the inflow conditions, for the v
     resolute can be specified as zero. However, this does modify
     the geometry near the inlet.
 
DO JJ=1,3
+ A=JJ-1;ZZ=ZWLAST*(A/NZ)**POWERZ
+ SETPT(1,1,JJ,0.0,0.0,ZZ);SETPT(2,1,JJ,1.0,0.0,ZZ)
ENDDO
 
   **For KK=4 on, the given analytical formula for the shape of
     the wall is employed
 
DO JJ=4,22
+ A=JJ-1;  ZZ=ZWLAST*(A/NZ)**POWERZ;  B=2-30*ZZ/RE
+ VV=2.7182818**B;  UU=1.0/VV;  DD=((VV-UU)/(VV+UU)-CT)/2.
+ SETPT(1,1,JJ,0.0,DD,ZZ)
+ SETPT(2,1,JJ,1.0,DD,ZZ)
ENDDO
   ** Set high boundary grid
DOMAIN(1,1,1,NY+1,NZ+1,NZ+1)
SETLIN(YC,YF+(YL-YF)*LNJ**POWERY)
   ** Set low boundary grid
DOMAIN(1,2,1,NY+1,1,1)
SETLIN(YC,YF+(YL-YF)*LNJ**POWERY)
   ** Set north boundary grid
DOMAIN(1,1,NY+1,NY+1,1,NZ+1)
SETLIN(ZC,ZF+(ZL-ZF)*LNK**POWERZ)
   ** Apply algebraic interpolation over domain
DOMAIN(1,2,1,NY+1,1,NZ+1);MAGIC(T)
   ** Apply Laplace solver for K greater than 2 to preserve
      shape of first two slabs for simple inlet conditions
SLIDN=T;JMON=10;KMON=10;MSWP=50
DOMAIN(1,2,1,NY+1,3,NZ+1);MAGIC(L);VIEW(I,1)
 
    GROUP 7.Variables stored, solved & named
SOLVE(P1,V1,W1);SOLUTN(P1,Y,Y,Y,N,N,N)
    GROUP 8. Terms (in differential equations) & devices
DIFCUT=0.0
    GROUP 9. Properties of the medium (or media)
ENUL=0.1
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=.5
    GROUP 13. Boundary conditions and special sources
   ** Inlet boundary conditions
DO II=1,21
+ Y0=((II-1)/NY)**POWERY;  Y1=(II/NY)**POWERY
+ YY=(Y1+Y0)/2;            C =3*(YY-YY**2/2)
+ INLET(INL:II:,LOW,1,1,II,II,1,1,1,1)
+ VALUE(INL:II:,P1,C*RHO1)
+ VALUE(INL:II:,W1,C)
ENDDO
   ** Outlet boundary conditions
PATCH(OUTLET,HIGH,1,1,1,NY,NZ,NZ,1,1);COVAL(OUTLET,P1,FIXP,0.071)
COVAL(OUTLET,V1,ONLYMS,0.0);COVAL(OUTLET,W1,ONLYMS,0.0)
   ** Friction at the South wall
WALL (SOUTH,SOUTH,1,1,1,1,1,NZ,1,1);COVAL(SOUTH,W1,1.0,0.0)
    GROUP 15. Termination of sweeps
LSWEEP=100
    GROUP 17. Under-relaxation devices
RELAX(V1,FALSDT,1.0); RELAX(W1,FALSDT,1.0)
    GROUP 22. Spot-value print-out
IZMON=10
    GROUP 23. Field print-out and plot control
ITABL=3;NPLT=2;NYPRIN=3;NZPRIN=3;TSTSWP=-1
SELREF=T; RESFAC=0.01
PATCH(INNER,PROFIL,1,1,1,1,1,NZ,1,1)
PLOT(INNER,P1,0.0,0.0);PLOT(INNER,W1,0.0,0.0)