PHOTON USE p phi msg EGWF TESTS IN 2/3D msg msg Pressure distribution: cl;red con p1 z max fi;.001 vec z max msg msg Pressto continue pause msg msg Temperature distribution: cl;red con tem1 z max fi;.001 vec z max msg msg Press to continue pause msg msg Concentration distribution: cl;red con c1 z max fi;.001 vec z max msg msg Press to continue pause msg msg Press e to END enduse GROUP 1. Run title and other preliminaries CHAR(CUSEGW,CTURB,CBFC);BOOLEAN(TURBL) INTEGER(NXF,NXL,NYF,NYL,NZF,NZL,ISUM,N1) REAL(VIN,VFRAC,TIN,CIN) VIN=10.;TIN=100.;CIN=-1. MESG( MESG( TEST OF EGWF FOR FLOW AROUND A HEATED SOLUBLE CUBE MESG( MESG( The flow enters the W-S-L corner of the domain, MESG( impinges on the solid cube, which is heated and MESG( which has constant C1 concentration, and then MESG( leaves the domain through the E-N-H corner. MESG( MESG( To run the EGWF variant, enter T MESG( To run the SOWF variant, enter F MESG( READVDU(CUSEGW,CHAR,F) MESG( MESG( Enter the turbulence model option; default MESG( a LAMinar flow. MESG( The options are: MESG( LAM - laminar flow MESG( TURB - simple mixing-length turbulence MESG( LVEL - LVEL algebraic turbulence model MESG( KLMODL - the k-l model of turbulence MESG( KEMODL - the k-e model of turbulence MESG( MESG( For further info. please see the PHOENICS MESG( encyclopaedia MESG( READVDU(CTURB,CHAR,LAM) MESG( MESG( Enter the COefficient for the WALL-functions; MESG( the default is GRND2 (LOGLAW). MESG( The options are: MESG( 1/PRNDTL() - suitable for laminar flow MESG( GRND1 (BLASIUS) - suitable for laminar flow MESG( GRND2 (LOGLAW) - equilibrium log-law wall- MESG( function for turbulent flows MESG( GRND3 (GENLAW) - non-equilibrium log-law wall- MESG( function for turbulent flows MESG( (not implemented for EGWF=T) MESG( READVDU(WALLCO,REAL,GRND2) MESG( MESG( Flow is always in the Y direction: MESG( by default the domain has 7x7x7 cells MESG( Enter the new NX, NY & NZ; MESG( if(iqalib.ne.0) then + nx=7; ny=7; nz=7 endif READVDU(NX,INT,7) READVDU(NY,INT,7) READVDU(NZ,INT,7) MESG( MESG( Is a BFC grid required? Enter T for True, MESG( N for False. The default is F. MESG( READVDU(CBFC,CHAR,F) MESG( MESG( Enter the frequency of block-corrections. MESG( The default is no block corrections. MESG( READVDU(ISOLBK,INT,0) EGWF=:CUSEGW:;BFC=:CBFC: IF (NX .LE. 0) THEN + NX=1 ENDIF IF (NY .LE. 0) THEN + NY=1 ENDIF IF (NZ .LE. 0) THEN + NZ=1 ENDIF IF (EGWF) THEN + CASE :CTURB: OF + WHEN TURB,4 + TEXT( TURBULENT TEST OF EGWF - EGWF VARIANT) + WHEN LAM,3 + TEXT( LAMINAR TEST OF EGWF - EGWF VARIANT) + ORELSE + TEXT( :CTURB: TEST OF EGWF - EGWF VARIANT) + TURBL = T + ENDCASE ELSE + CASE :CTURB: OF + WHEN TURB,4 + TEXT( TURBULENT TEST OF EGWF - SOWF VARIANT) + WHEN LAM,3 + TEXT( LAMINAR TEST OF EGWF - SOWF VARIANT) + ORELSE + TEXT( :CTURB: TEST OF EGWF - SOWF VARIANT) + TURBL = T + ENDCASE ENDIF GROUP 3. X-direction grid specification GRDPWR(X,NX,1.0,1.0) NXF=NX/3+1;NXL=NX-NX/3 GROUP 4. Y-direction grid specification GRDPWR(Y,NY,1.0,1.0) NYF=NY/3+1;NYL=NY-NY/3 GROUP 5. Z-direction grid specification GRDPWR(Z,NZ,1.0,1.0) NZF=NZ/3+1;NZL=NZ-NZ/3 GROUP 7. Variables stored, solved & named SOLVE(P1);SOLUTN(P1,Y,Y,Y,P,P,Y) SOLVE(TEM1);SOLUTN(TEM1,Y,Y,Y,P,P,Y) SOLUTN(C1,Y,Y,Y,P,P,Y) ISUM=0 IF (NX .GT. 1) THEN + SOLVE(U1);ISUM=ISUM+1 + SOLUTN(U1,P,P,Y,P,P,Y) ENDIF IF (NY .GT. 1) THEN + SOLVE(V1);ISUM=ISUM+1 + SOLUTN(V1,P,P,Y,P,P,Y) ENDIF IF (NZ .GT. 1) THEN + SOLVE(W1);ISUM=ISUM+1 + SOLUTN(W1,P,P,Y,P,P,Y) ENDIF VFRAC=ISUM**0.5 STORE(PRPS) IF (ISOLBK .NE. 0) THEN + STORE(BLOK) + IVARBK=-1 ENDIF IF (TURBL) THEN + TURMOD(:CTURB:) + SOLUTN(KE,P,P,p,P,P,Y) + SOLUTN(EP,P,P,p,P,P,Y) ELSE IF (:CTURB: .EQ. TURB) THEN + STORE(LEN1,VIST) ENDIF IF (EGWF) THEN + STORE(SKIN,STAN,STRS) ENDIF GROUP 8. Terms (in differential equations) & devices *** Get rid of potentially confusing heat sources that are *** irrelevant to the tests being performed TERMS(TEM1,N,P,P,P,P,P) *** Use the upwind scheme to get more easily understandable *** behaviour DIFCUT=0.0 GROUP 9. Properties of the medium (or media) *** Set up the properties of the materials to be used *** NB. the viscosity is boosted to make the viscous *** effects more noticeable CSG10='q1' MATFLG=T;NMAT=2 55 1.0 1.e-3 1000. 2.5 4.0e-3 155 8000. 1.0 500. 50. 0.0 *** Set up the turbulence models for the turbulent cases *** NB. at present limited to an artificially simple mixing *** length model IF (:CTURB: .EQ. TURB .OR. :CTURB: .EQ. KLMODL) THEN + IF (NX.GT.1) THEN + EL1=LINEARX + ELSE + EL1=LINEARY + ENDIF + EL1A=(XULAST**2+YVLAST**2+ZWLAST**2)**0.5;EL1B=0.0 ENDIF IF (:CTURB: .EQ. TURB) THEN + ENUT=PROPLEN;ENUTA=0.0;ENUTB=0.1 ENDIF GROUP 11. Initialization of variable or porosity fields isum vfrac FIINIT(U1)=VIN/VFRAC;FIINIT(V1)=VIN/VFRAC;FIINIT(W1)=VIN/VFRAC FIINIT(PRPS)=55;PRNDTL(TEM1)=CONDFILE;ENUL=FILE; FIINIT(TEM1)=10.;FIINIT(C1)=0.0 IF (.NOT. EGWF) THEN + CONPOR(BLOCKCP,-1,CELL,-NXF,-NXL,-NYF,-NYL,-NZF,-NZL) ENDIF CASE :CTURB: OF WHEN KEMODL,6 + FIINIT(KE) = 0.01*VIN*VIN + FIINIT(EP) = FIINIT(KE)**1.5 WHEN KLMODL,6 + FIINIT(KE) = 0.01*VIN*VIN ENDCASE INIADD=F PATCH(BLOCK ,CELL,NXF,NXL,NYF,NYL,NZF,NZL,1,1) PATCH(BLOCKI,INIVAL,NXF,NXL,NYF,NYL,NZF,NZL,1,1) INIT (BLOCKI,PRPS,0.0,155.) IF (ISOLBK .NE. 0) THEN + FIINIT(BLOK)=1 + INIT(BLOCKI,BLOK,0.0,2.0) ENDIF IF (ISOLBK .NE. 0) THEN + FIINIT(BLOK)=1.0;INIT(BLOCKI,BLOK,0.0,2.0) ENDIF IF (EGWF) THEN + FIINIT(SKIN)=0.0;FIINIT(STAN)=0.0 ENDIF GROUP 13. Boundary conditions and special sources *** Set the C1 and TEM1 sbources in the blocks: fixed *** flux for heat(??), and fixed values of C1 (heated *** rock-salt dissolving into water??) COVAL(BLOCK,C1,FIXVAL,-CIN) COVAL(BLOCK,TEM1,FIXFLU,100.) *** Setup inlet conditions IF (NX .GT. 1) THEN + INLET(XINL,WEST,1,1,1,NYF-1,1,NZF-1,1,LSTEP) + VALUE(XINL,P1,RHO1*VIN) ;VALUE(XINL,U1,VIN/VFRAC) + VALUE(XINL,V1,VIN/VFRAC);VALUE(XINL,W1,VIN/VFRAC) + VALUE(XINL,TEM1,TIN);VALUE(XINL,C1,CIN) + OUTLET(XOUTL,EAST,NX,NX,NYL+1,NY,NZL+1,NZ,1,LSTEP) + VALUE(XOUTL,P1,0.0) ENDIF IF (NZ .GT. 1) THEN + INLET(ZINL,LOW,1,NXF-1,1,NYF-1,1,1,1,LSTEP) + VALUE(ZINL,P1,RHO1*VIN) ;VALUE(ZINL,U1,VIN/VFRAC) + VALUE(ZINL,V1,VIN/VFRAC);VALUE(ZINL,W1,VIN/VFRAC) + VALUE(ZINL,TEM1,TIN);VALUE(ZINL,C1,CIN) + OUTLET(ZOUTL,HIGH,NXL+1,NX,NYL+1,NY,NZ,NZ,1,LSTEP) + VALUE(ZOUTL,P1,0.0) ENDIF IF (NY .GT. 1) THEN + INLET(YINL,SOUTH,1,NXF-1,1,1,1,NZF-1,1,LSTEP) + VALUE(YINL,P1,RHO1*VIN) ;VALUE(YINL,U1,VIN/VFRAC) + VALUE(YINL,V1,VIN/VFRAC);VALUE(YINL,W1,VIN/VFRAC) + VALUE(YINL,TEM1,TIN);VALUE(YINL,C1,CIN) + OUTLET(YOUTL,NORTH,NXL+1,NX,NY,NY,NZL+1,NZ,1,LSTEP) + VALUE(YOUTL,P1,0.0) ENDIF GROUP 15. Termination of sweeps LSWEEP=50 GROUP 16. Termination of iterations LITER(TEM1)=20;LITER(C1)=20;LITER(U1)=10;LITER(V1)=10 GROUP 17. Under-relaxation devices RELAX(U1,FALSDT,5.*XULAST/(NX*VIN)) RELAX(V1,FALSDT,5.*YVLAST/(NY*VIN)) RELAX(W1,FALSDT,5.*ZWLAST/(NZ*VIN)) RELAX(KE,FALSDT,5.*YVLAST/(NY*VIN)) RELAX(EP,FALSDT,5.*YVLAST/(NY*VIN)) RELAX(TEM1,FALSDT,1000.) RELAX(C1, FALSDT,1000.) GROUP 18. Limits on variables or increments to them IF (EGWF) THEN + VARMIN(SKIN)=0.0;VARMIN(STAN)=0.0 ENDIF GROUP 21. Print-out of variables OUTPUT(P1,Y,P,P,Y,Y,Y);OUTPUT(U1,Y,P,P,Y,Y,Y) OUTPUT(V1,Y,P,P,Y,Y,Y);OUTPUT(W1,Y,P,P,Y,Y,Y) OUTPUT(TEM1,Y,P,P,Y,Y,Y);OUTPUT(C1,Y,P,P,Y,Y,Y) IF (:CTURB: .EQ. KEMODL .OR. :CTURB: .EQ. KLMODL) THEN + OUTPUT(KE,Y,P,P,Y,Y,Y) ENDIF IF (BFC) THEN + IF (NX .GT. 1) THEN + OUTPUT(UCRT,Y,P,P,Y,Y,Y) + ENDIF + IF (NY .GT. 1) THEN + OUTPUT(VCRT,Y,P,P,Y,Y,Y) + ENDIF + IF (NZ .GT. 1) THEN + OUTPUT(WCRT,Y,P,P,Y,Y,Y) + ENDIF ENDIF GROUP 22. Spot-value print-out IXMON=NXF;IYMON=NY/2+1;IZMON=NZ/2+1 TSTSWP=-1 GROUP 23. Field print-out and plot control NXPRIN=1;NYPRIN=1;NZPRIN=1 IF (NX .GT. 1 .AND. NY .GT. 1) THEN + PATCH(CONZ,CONTUR,1,NX,1,NY,NZ/2+1,NZ/2+1,1,1) + PLOT (CONZ,P1,0.0,20.) + PLOT (CONZ,TEM1,0.0,20.) + PLOT (CONZ,C1,0.0,20.) ENDIF IF (NY .GT. 1 .AND. NZ .GT. 1) THEN + PATCH(CONX,CONTUR,NX/2+1,NX/2+1,1,NY,1,NZ,1,1) + PLOT (CONX,P1,0.0,20.) + PLOT (CONX,TEM1,0.0,20.) + PLOT (CONX,C1,0.0,20.) ENDIF IF (NZ .GT. 1 .AND. NX.GT. 1) THEN + PATCH(CONY,CONTUR,1,NX,NY/2+1,NY/2+1,1,NZ,1,1) + PLOT (CONY,P1,0.0,20.) + PLOT (CONY,TEM1,0.0,20.) + PLOT (CONY,C1,0.0,20.) ENDIF IF (NX .GT. 1) THEN + PATCH(PRFX,PROFIL,1,NX,NY/2+1,NY/2+1,NZ/2+1,NZ/2+1,1,1) + PLOT (PRFX,TEM1,100.0,1000.) + PLOT (PRFX,C1,-1.,1.) ENDIF IF (NY .GT. 1) THEN + PATCH(PRFY,PROFIL,NX/2+1,NX/2+1,1,NY,NZ/2+1,NZ/2+1,1,1) + PLOT (PRFY,TEM1,100.0,1000.) + PLOT (PRFY,C1,-1.,1.) ENDIF IF (NZ .GT. 1) THEN + PATCH(PRFZ,PROFIL,NX/2+1,NX/2+1,NY/2+1,NY/2+1,1,NZ,1,1) + PLOT (PRFZ,TEM1,100.0,1000.) + PLOT (PRFZ,C1,-1.,1.) ENDIF