```
PHOTON USE
p

gr x 1
msg Grid
vec x 1 sh
MSG Velocity vectors
pause
gr cl; gr ou x 1; red
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)
```