Author's explanation 01.04.09 ---------------------------------------- This case is a development of case 290 which concerns turbulent flow over and downstream of a backward-facing step but uses a grid which is far too coarse for accuracy. It has therefore been modified so as to enable the effects of the grid-defining parameters NX, NY, XPOWER and YPOWER to be easily varied. This is done by changing the name of the STEP Patch to .STEP . The '.' which starts its name is a signal to EARTH that locating arguments are geometric rather than indicial, albeit non- dimensionalised by reference to the domain dimensions. [The following statement was true when the .patch was first introduced, but it the multiply-by-1000 practice is no longer needed: A multiplier of 1000 is used (and matched by division on arrival in EARTH) because, although the arguments are real numbers, the SATELLITE reduces them to whole numbers.] Although it is not necessary, because they lie at domain boundaries, the inlet, outlet and wall boundary conditions are also expressed via dot-patches so as to test and illustrate the technique. DBS ----------------------------------------------------------------- Setting the grid fineness and distribution interactively NX=20; ny=15 ! defaults real(xpower,ypower) xpower=1.0; ypower=1.0 if(iqalib.ne.0) then + ans=y endif label message mesgm(nx=:nx: ny=:ny: xpower=:xpower: ypower=:ypower: If ok, RETURN readvdu(ans,char,y) if(:ans:.eq.y) then goto begin else mesg(insert required NX, or RETURN readvdu(nx,int,nx) mesg(insert required NY, or RETURN readvdu(ny,int,ny) mesg(insert required XPOWER, or RETURN readvdu(xpower,real,xpower) mesg(insert required YPOWER, or RETURN readvdu(ypower,real,ypower) goto message endif label begin 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 REAL(HEIGHT,WIDTH,CLEN,SLEN,UIN,TKEIN,EPSIN,GENUT) UIN=1.0 ** Calculation of domain specifications HEIGHT=0.0381; WIDTH=3.*HEIGHT SLEN=4.*HEIGHT; CLEN=20.*HEIGHT XULAST=CLEN YVLAST=WIDTH GRDPWR(X,NX,XULAST,XPOWER) GRDPWR(Y,NY,YVLAST,YPOWER) 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 ep0 ENDIF GROUP 9. Properties of the medium (or media) 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 EGWF=T;STORE(PRPS) declarations and settings of dot-patch arguments REAL(RXF,RXL,RYF,RYL,RZF,RZL) The next three de-activated statements illustrate the now discarded multiplication-by-1000 practice. RXF=0.0; RXL=SLEN/XULAST*1000 RYF=0.0; RYL=HEIGHT/YVLAST*1000 RZF=0.0; RZL=ZWLAST/ZWLAST*1000 RXF=0.0; RXL=SLEN/XULAST RYF=0.0; RYL=HEIGHT/YVLAST RZF=0.0; RZL=ZWLAST/ZWLAST PATCH(.STEP,INIVAL,RXF,RXL,RYF,RYL,RZF,RZL,1,1) ! Important warning ! Although it might seem reasonable to do so, DO NOT set the ! last but one arguments above to 0 because doing so confuses the ! indexing system of the satellite, with unforeseeable results. ! Note also that 1,1000 is satisfactory for the last two ! arguments of the group 13 dot-patches; but the equally ! reasonable 0,1000, although OK for the x, y and z arguments, ! causes the aforesaid damage. COVAL(.STEP,PRPS,0.0,198.0) TEXT( :char1: model; EGWF= :EGWF: TITLE GROUP 13. Boundary conditions and special sources ** Inlet PATCH(INLET,WEST,1,1,1,NY,1,1,1,1) ! the no-dot patch PATCH(.INLET,WEST,0,0,0,1000,0,1000,1,1000) ! the discarded ! multiply-by-1000 practice PATCH(.INLET,WEST,0,0,0,1,0,0,0,1) COVAL(.INLET,P1,FIXFLU,RHO1*UIN) COVAL(.INLET,U1,ONLYMS,UIN) COVAL(.INLET,KE,ONLYMS,TKEIN) COVAL(.INLET,EP,ONLYMS,EPSIN) ** Exit PATCH(OUTLET,EAST,NX,NX,1,NY,1,1,1,1) ! the no-dot patch PATCH(.OUTLET,EAST,1000,1000,0,1000,0,1000,1,1000) ! the discarded ! multiply-by-1000 practice PATCH(.OUTLET,EAST,1,1,0,1,0,1,0,1) COVAL(.OUTLET,P1,1.0E5,0.0) ** N-wall WALL (WFNN,NORTH,1,NX,NY,NY,1,1,1,1) ! the no-dot patch WALL (.WFNN,NORTH,0,1000,1000,1000,0,1000,1,1000) ! the discarded ! multiply-by-1000 practice WALL (.WFNN,NORTH,0,1,1,1,0,1,1,1) ** S-wall WALL (.WFNS,SOUTH,0,1000,0,0,0,1000,1,1000) ! the discarded ! multiply-by-1000 practice WALL (.WFNS,SOUTH,0,1,0,0,0,1,0,1) GROUP 15. Termination of sweeps LSWEEP=1000;SELREF=T;RESFAC=0.01 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=NY/2;IXMON=NX/2;NPRMON=100;TSTSWP=12;UWATCH=T 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) tstswp=-1 lsweep=1 debug=t dbindx=t search=t idbf0=414 flag=t stop