TALK=T;RUN(1,1)
dbs 21.05.10 provision made (ixfac,iyfac) for increasing
grid fineness
290
PHOTON USE
p;;;;
use patgeo
gr ou z 1
msg contours of velocity, and velocity vectors
con u1 z 1 fi;0.01;vec z 1
msg press RETURN for length scale
pause;con off;vec off;red
msg contours of turbulence length scale
con len1 z 1 fi;0.01;pause;con off;red
msg contours of effective viscosity
con enut z 1 fi;0.01
enduse
GROUP 1. Run title and other preliminaries
TEXT(Backward-Facing Step; K-E Or LVEL
TITLE
DISPLAY
This case involves steady, turbulent flow over a backward-facing
step, of height h, in a two-dimensional channel of width 3*h.
The calculations are started at a distance 4h upstream of the step
and terminated at a distance 16h downstream.
//////////////////////// wall ///////////////////////
-----------------------------------------------------
Pressure
Inlet -------> ---------> fixed at
zero
/________________
| ////////////////| Exit
| wall /| Recirculation
| /| <---- ----->
\ constant /|____________________________________
prescribed mass //////////////////////////////////////
inflow rate wall
ENDDIS
INTEGER(NYS,NXS)
REAL(HEIGHT,WIDTH,CLEN,SLEN,REYNO,UIN,TKEIN,EPSIN,GENUT)
** Calculation of domain specifications
HEIGHT=0.0381;WIDTH=3.*HEIGHT;SLEN=4.*HEIGHT;CLEN=20.*HEIGHT
NXS=4;NYS=10;REYNO=4.5E4;UIN=13.0
mesg(Reynolds number based on step height and inlet velocity
mesg(is :reyno: OK? If not, insert other value
readvdu(reyno,real,reyno)
integer(ixfac,iyfac) ! increase ixfac and iyfac in order to
! use finer grid
ixfac=4;iyfac=4
ixfac=1;iyfac=1
nxs=nxs*ixfac
nys=nys*iyfac
GROUP 3. X-direction grid specification
** Full length of channel = 0.762
NREGX=3
IREGX=1;GRDPWR(X,NXS,SLEN,1.0)
IREGX=2;GRDPWR(X,8*ixfac,0.25*(CLEN-SLEN),1.0)
IREGX=3;GRDPWR(X,8*ixfac,0.75*(CLEN-SLEN),1.0)
GROUP 4. Y-direction grid specification
** Full width of channel = 0.1143
NREGY=2
IREGY=1;GRDPWR(Y,NYS,HEIGHT,1.0)
IREGY=2;GRDPWR(Y,5*iyfac,WIDTH-HEIGHT,1.0)
GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1);SOLUTN(P1,Y,Y,Y,N,N,N);STORE(ENUT,LEN1)
mesg(Turbulence model is KEMODL. Use LVEL instead? (y/n)
CHAR(char1)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
TURMOD(LVEL)
char1=LVEL
store(wgap)
ELSE
TURMOD(KEMODL)
char1=KE-EPS
KELIN=3
mesgb(kelin=3; no under-relaxation is used for ke or ep
ENDIF
GROUP 9. Properties of the medium (or media)
RHO1=1.0;ENUL=UIN*HEIGHT/REYNO
GROUP 11. Initialization of variable or porosity fields
** Calculation of KE (where fric=0.018)...
TKEIN=0.25*UIN*UIN*0.018
** Calculation of EP (where lmix=0.09 x h)...
EPSIN=TKEIN**1.5*0.1643/(0.09*height)
** Initial values
FIINIT(U1)=UIN;FIINIT(P1)=1.3E-4;FIINIT(KE)=TKEIN;FIINIT(EP)=EPSIN
** Initialization of variables in blocked region
mesga(Use EARTH-generated wall functions? (y/n)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
EGWF=T;STORE(PRPS)
PATCH(STEP,INIVAL,#1,#1,#1,#1,#1,#1,1,1)
COVAL(STEP,PRPS,0.0,198.0)
ELSE
EGWF=F;CONPOR(STEP,0.0,CELL,#1,-#1,#1,-#1,#1,#1)
ENDIF
TEXT( :char1: model; Re = :reyno:; EGWF= :EGWF:
TITLE
GROUP 13. Boundary conditions and special sources
** Inlet
INLET(INLET,WEST,#1,#1,#2,#NREGY,#1,#1,1,1)
VALUE(INLET,P1,UIN);VALUE(INLET,U1,UIN)
VALUE(INLET,KE,TKEIN);VALUE(INLET,EP,EPSIN)
** Exit
PATCH(OUTLET,EAST,#NREGX,#NREGX,#1,#NREGY,#1,#1,1,1)
COVAL(OUTLET,P1,1.0E5,0.0)
COVAL(OUTLET,U1,ONLYMS,0.0);COVAL(OUTLET,V1,ONLYMS,0.0)
COVAL(OUTLET,KE,ONLYMS,0.0);COVAL(OUTLET,EP,ONLYMS,0.0)
** N-wall
WALL (WFNN,NORTH,#1,#NREGX,#NREGY,#NREGY,#1,#1,1,1)
** S-wall
WALL (WFNS,SOUTH,#2,#NREGX,#1,#1,#1,#1,1,1)
GROUP 15. Termination of sweeps
LSWEEP=100;SELREF=T;RESFAC=0.01
lsweep=lsweep*2
isg21=lsweep ! to endure that the run does not terminate
! prematurely
GROUP 16. Termination of iterations
LITER(U1)=20;LITER(V1)=20
GROUP 17. Under-relaxation devices
RELAX(U1,FALSDT,SLEN/UIN);RELAX(V1,FALSDT,SLEN/UIN)
relax(ke,linrlx,0.5);relax(ep,linrlx,0.5)
GROUP 22. Monitor print-out
IYMON=NYS-2;IXMON=NXS+2;NPRMON=100;TSTSWP=12;UWATCH=T
tstswp=-1 ! to activate the graphical monitor
#maxmin ! to cause maximum and minimum values to be
! plotted on graphical monitor
GROUP 23. Field print-out and plot control
NPLT=1;IPLTL=LSWEEP;ITABL=1
PATCH(MAP,CONTUR,1,NX,1,NY,1,1,1,1)
PLOT(MAP,P1,0.0,10.0);PLOT(MAP,U1,0.0,10.0);PLOT(MAP,V1,0.0,10.0)