** LOAD(105) from the PHOENICS Input Library GROUP 1. Run title and other preliminaries TEXT(1D TRANS FOR H1 WITH FLOW + SOURCE #cls TITLE DISPLAY This problem deserves study; for it exemplifies, in a simple form, the main ingredients of computational fluid dynamics, namely the effects of: * time-dependence, * convective transport * diffusion (laminar and turbulent) * sources, and * non-uniform properties. The problem-defining data are set in two alternative ways, namely: 1. using INIT and COVAL for initial and boundary conditions, and GRND settings for the properties; and 2. using In-Form statements, which are easier to understand, and provide to freedom to set whatever formulae are desired. The results are of course the same, as are (near enough) the computer times. PHOTON USE commands are provided for ease of display. These exploit the possibility of treating the time dimension as though it were z-direction distance. It is the command IDISPA=1 which effects this. ENDDIS boolean(inform) mesg(In-Form for setting data ? (Y/n) readvdu(ans,char,y) inform = :ans:.eq.y inform GROUP 2. Transience; time-step specification STEADY=F;GRDPWR(T,20,1.0E-4,1.0) GROUP 3. X-direction grid specification GRDPWR(X,11,0.11,1.0) GROUP 7. Variables stored, solved & named Storage is provided for both laminar and turbulent viscosities STORE(ENUL,ENUT,TMP1,EL1) **First-phase enthalpy is the variable solved . SOLVE(H1) **The first-phase velocity is allocated storage, because it is desired that convective effects shall be present, but U1 is not solved. STORE(U1) GROUP 8. Terms (in differential equations) & devices **The built-in source for enthalpy is switched off by the following N. It is the Dp/Dt term in the regular enthalpy equation, which is not appropriate because pressure is not being solved for in this calculation. TERMS(H1,N,Y,Y,Y,Y,Y) GROUP 9. Properties of the medium (or media) **The following statements make temperature a linear function of enthalpy (H1), indeed equal to it. if(inform) then (stored var tmp1 is h1) else TMP1=LINH; TMP1A=0.0; CP1=1.0 endif **The following statements make the length scale a "ramp" function of distance x. See GREX3 Group 9 Section 12. if(inform) then (property var el1 is min((0.0 + 1.0*xg) , 0.05)) else EL1=LINEARX; EL1A=0.0; EL1B=1.0; EL1C=0.05 endif **The following statements make the laminar viscosity equal 0.5+0.5*temperature . if(inform) then (property var enul is 0.5 + 0.5*tmp1) else ENUL=LINTEM; ENULA=0.5; ENULB=0.5 endif **The following statements make ENUT equal 0.0+100.0*length-scale if(inform) then (property var enut is 0.0 + 100.0 * el1) else ENUT=PROPLEN;ENUTA=0.0;ENUTB=100.0 endif GROUP 11. Initialization of variable or porosity fields ** Of the following initial values, only that of U1 will prevail throughout the computation. FIINIT(H1)=0.5; FIINIT(U1)=1.E2; FIINIT(ENUL)=1.0; FIINIT(ENUT)=1.0 FIINIT(TMP1)= 0.5 ** Introduce a step into the initial enthalpy profile, which will otherwise have the value FIINIT(H1) everywhere, over the IX range from 3 to 8 inclusive. PATCH(INI1,INIVAL,3,8,1,1,1,1,1,1) IF(INFORM) THEN (INITIAL OF H1 AT INI1 IS 1.0) ELSE INIT(INI1,H1,0.0,1.0) ENDIF GROUP 13. Boundary conditions and special sources **Enthalpy is held to zero at the low-x end, PATCH(X1,CELL,1,1,1,NY,1,NZ,1,LSTEP) if(inform) then (SOURCE of H1 at X1 is 0.0 with FIXV) else COVAL(X1,H1,FIXVAL,0.0) endif **Enthalpy is held to 1.0 at the high-x end PATCH(XNX,CELL,NX,NX,1,NY,1,NZ,1,LSTEP); if(inform) then (SOURCE of H1 at XNX is 1.0 with FIXV) else COVAL(XNX,H1,FIXVAL,1.0) endif ** A volumetric enthalpy source is supplied over the whole of the integration domain PATCH(WHOLE,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP) if(inform) then (SOURCE of H1 at WHOLE is 1.e4*(2.0-h1) with LINE) else COVAL(WHOLE,H1,1.E4,2.0) endif GROUP 15. Termination of sweeps **Because of the special nature of the problem specification in which convection is being introduced within the domain, but there is no mass flow communication with outside, it can be expected that large "residuals" will be reported. However, RESREF must be set; otherwise iteration will terminate prematurely. LSWEEP will control the termination in this case. LSWEEP=100;RESREF(H1)=1.0 SPEDAT(SET,GXMONI,TRANSIENT,L,F) GROUP 21. Print-out of variables **Provide printout of the two components of the reference viscosity OUTPUT(ENUL,Y,N,N,N,Y,Y); OUTPUT(ENUT,Y,N,N,N,Y,Y) GROUP 23. Field print-out and plot control **Tabulate spot-values and residuals ITABL=2 **Print every second value, starting at IX=2 and ending at IX=10. NXPRIN=2; IXPRF=2; IXPRL=10 **Plot profiles over the length of the domain PATCH(XWISE,PROFIL,1,NX,1,1,1,1,1,LSTEP) PLOT(XWISE,H1,0.0,0.0); PLOT(XWISE,ENUL,0.0,0.0) PLOT(XWISE,ENUT,0.0,0.0) **Plot a profile of the values at IX=5 with time as the abscissa PATCH(TIMEWISE,PROFIL,5,5,1,1,1,1,1,LSTEP) PLOT(TIMEWISE,H1,0.0,0.0); PLOT(TIMEWISE,ENUL,0.0,0.0) LIBREF = 105 TSTSWP=-1 IDISPA=1 PHOTON USE p parphi 1 1 1.e3 gr ou y 1 msg(vertical coord. = distance, x; horizontal coord. = time pause con h1 y 1 fi;0.001 pause con enut y 1 fi;0.001 pause con enul y 1 fi;0.001 ENDUSE