GROUP 1. Run title and other preliminaries
TEXT(YX LAMINAR BACKWARD-FACING-STEP FLOW
TITLE
mesg( PC486/50 time last reported as 6.0 min
DISPLAY
Numerical Schemes validation example:
2-d x-y, Cartesian, steady, elliptic simulation
The case considered is steady incompressible, laminar backward-
facing-step flow, i.e. flow through a straight channel having a
sudden asymmetric expansion. The flow is characterised by the
step height H and the Reynolds number Re based on the bulk inlet
inlet velocity and 2h, where h is the flow inlet height. The
channel expansion ratio is 1.94, and the total length of the
domain is 40H. A fully-developed parabolic velocity profile is
prescribed at the inlet boundary located at the step. The
calculation is made for Re=150 with the cubic upwind scheme, and
the option exists to select the hybrid or van Albada scheme.
ENDDIS
This test problem is widely used for assessing the accuracy of
numerical methods because of the dependence of the reattachment
lengths X on Re. The flow has been studied experimentally by
Armaly et al (J.Fluid Mech., Vol.127, p473 (1983)), and
numerically in several papers, e.g.Gartling (Int.J.Num.Meth.Fluids
Vol.11, p953 (1990)) and Freitas(ASME J.Fluids Engng, Vol.117,
p208, (1995)), Gresho et al (Int.J.Num.Meth.Fluids Vol.17, p501
(1993))).
The user is invited to perform 2 calculations for comparison with
Armaly's data, namely Re=150 and Re=450. At Re=150, a primary
recirculation zone develops immediately downstream of the step
with X1/H=5.0. At Re=450, X1/H=9.5 and an additional separation
cell forms on the upper wall of the channel. The measurements
indicate that the separation point of this cell is located at
X2/H=7.6 and the reattachment point at X3/H=11.3, yielding a cell
length of DX/H=3.7. For Re > 450, the flow begins to exhibit 3d
effects. The main results may be summarised as follows:
Re=150 Data Hybrid Cubic-Upwind Van-Albada
X1/H 4.2 4.17 4.24 4.25
Re=450
X1/H 9.5 8.79 9.13 9.06
X2/H 7.6 7.66 8.26 8.09
X3/H 11.3 11.14 11.39 11.52
DX/H 3.7 3.48 3.13 3.43
Although no grid-refinement studies have been performed, the
results are in fairly good agreement with the data. It should
also be mentioned that the foregoing results are better than
those reported by Freitas(1995) for other general-purpose codes.
The measurements indicate that for Re>450, 3d effects become
significant, and for Re >6,600 the flow is 2d fully-turbulent.
PHOTON USE
P
0.20443E+04 0.15633E+04 CR
gr ou z 1;mag gr 5
0.28838E+03 0.17522E+04 CR
msg Reynolds number = 150 Cubic upwind scheme
vec z 1 x 1 40 y 1 m sh
STREAM 2D Z 1 X 1 40 Y 1 M
-.699 0. 5
STREAM 2D Z 1 X 1 40 Y 1 M
0. 12. 8
msg press to continue
msg press to end
pause
ENDUSE
AUTOPLOT USE
file
phi 5
d 1 u1 y 1;plot;div x .49 1
scale x 0 6;level y 0
msg Reynolds number = 150 Cubic upwind scheme
msg horizontal (U1) velocity distribution along bottom wall
msg reattachment point is where U1 crosses zero point
msg press to continue
msg press to end
pause
ENDUSE
CHAR(SCHM);REAL(UIN,DY,HIN,HSTEP,YH,LENGTH,LENX1,RENO,DUM)
REAL(UINAV,RLXFAC);INTEGER(NX1,NX2,NY1,NY2,NY2P1)
** All dimensions are based on: g, cm, sec
HIN=0.52;HSTEP=0.49;LENGTH=40.0*HSTEP
MESG(Enter Reynolds number - default 150
READVDU(RENO,REAL,150.)
ENUL=0.16;UINAV=RENO*ENUL/(2.0*HIN)
UINAV
GROUP 2. Transience; time-step specification
GROUP 3. X-direction grid specification
was nx1=200 nx2=50
NX1=160;NX2=40;LENX1=0.666666*LENGTH
NREGX=2;IREGX=1;GRDPWR(X,NX1,LENX1,1.0)
IREGX=2;GRDPWR(X,NX2,LENGTH-LENX1,1.2)
GROUP 4. Y-direction grid specification
NY2=16;NY1=16;NREGY=2;IREGY=1;GRDPWR(Y,NY1,HSTEP,1.0)
IREGY=2;GRDPWR(Y,NY2,HIN,1.0);DY=HIN/NY2
GROUP 5. Z-direction grid specification
NZ=1;ZWLAST=0.01
GROUP 6. Body-fitted coordinates or grid distortion
GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1)
GROUP 8. Terms (in differential equations) & devices
MESG( Enter required convection scheme
MESG( Default: CUS - Cubic upwind
MESG( The options are:
MESG( CUS - Cubic upwind
MESG( HYB - Hybrid scheme
MESG( VAN - VAN ALBADA scheme
READVDU(SCHM,CHAR,CUS)
CASE :SCHM: OF
WHEN CUS,3
+ MESG(Cubic upwind scheme
+ SCHEME(CUS,U1,V1)
WHEN FOU,3
+ MESG(First-order upwind scheme
+ DIFCUT=0.5
WHEN VAN,3
+ MESG(Van Albada scheme
+ SCHEME(VANALB,U1,V1)
ENDCASE
GROUP 9. Properties of the medium (or media)
RHO1=1.2E-3
GROUP 10. Inter-phase-transfer processes and properties
GROUP 11. Initialization of variable or porosity fields
FIINIT(U1)=0.5*UINAV
GROUP 12. Convection and diffusion adjustments
GROUP 13. Boundary conditions and special sources
** Inlet velocity profile: u(y)/uav = {6y(hin-y)}/hin**2
umax/uav=1.5 at y=hin/2. The profile presumes a uniform
mesh distribution.
REAL(SAIN,SFIN,UBAR,FCON);FCON=RHO1*DY*ZWLAST;SAIN=0.;SFIN=0.
DO JJ=1,NY2
+ YH=0.5*DY+DY*(JJ-1);UIN=6.0*UINAV*YH*(HIN-YH)
+ UIN=UIN/(HIN*HIN);SAIN=SAIN+FCON;SFIN=SFIN+UIN*FCON
+ PATCH(IN:JJ:,WEST,1,1,NY1+JJ,NY1+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
** check on ubar
UBAR=SFIN/SAIN
** Wall boundary conditions
PATCH(TWALL,NWALL,1,NX,NY,NY,1,NZ,1,LSTEP);COVAL(TWALL,U1,1.0,0.0)
PATCH(BWALL,SWALL,1,NX,1,1,1,NZ,1,LSTEP);COVAL(BWALL,U1,1.0,0.0)
PATCH(STPWL,WWALL,1,1,1,NY1,1,NZ,1,LSTEP);COVAL(STPWL,V1,1.0,0.0)
** Outlet boundary
PATCH(OUT,EAST,NX,NX,1,NY,1,NZ,1,LSTEP);COVAL(OUT,P1,10.0,0.0)
GROUP 14. Downstream pressure for PARAB=.TRUE.
GROUP 15. Termination of sweeps
** The number of sweeps required for convergence depends on the
scheme used. Typically, 1000-2000 sweeps are required with
the current mesh density.
LSWEEP=100
GROUP 16. Termination of iterations
LITER(U1)=10;LITER(V1)=10
GROUP 17. Under-relaxation devices
** The cubic-upwind & van-Albada schemes require RLXFAC
to be halved as convergence is approached.
RLXFAC=XULAST/(NX*UINAV)
RELAX(P1,LINRLX,1.0);RELAX(U1,FALSDT,1.0*RLXFAC)
RELAX(V1,FALSDT,1.0*RLXFAC)
GROUP 18. Limits on variables or increments to them
GROUP 22. Spot-value print-out
TSTSWP =-1;IYMON=NY-3;IXMON=130
GROUP 23. Field print-out and plot control
NPLT=20;ITABL=3
GROUP 24. Dumps for restarts