```
TEXT(Decay Of Turbulence. K-Epsilon Model
TITLE
DISPLAY
Homogeneous turbulence decays behind a grid in a stream of uniform
velocity. The transport equations of the k-epsilon model reduce to
ordinary differential equations with analytical solutions of which
that for k is:

k=( k0**a - a*e0*(x-x0)/{u0*k0**(1-a)} ) ** (1/a)

where k0 and e0 are k and e for x=0, and u0 is the flow velocity;
a=1.0-C2E. The value of C2E that fits the experimental data is 1.9

The analytical solution is computed in the Q1 file, by way of a DO
loop, and stored in the variable ANKE. The small differences
between the computed KE and ANKE result from the upwind
differencing: ANKE is the cell-centre value, whereas KE is the
convected-out value.
ENDDIS
REAL(XLEN,UIN,REYNO,WIDTH,TKEIN,EPIN)
XLEN=1.0;REYNO=5.E5;UIN=10.;WIDTH=0.25
GROUP 3. X-direction grid specification
GRDPWR(X,200,XLEN,1.0)
GROUP 7. Variables stored, solved & named
STORE(U1,ENUT);TURMOD(KEMODL);NAME(ENUT)=VIST
GROUP 8. Terms (in differential equations) & devices
DIFCUT=0.0;ENUL=UIN*WIDTH/REYNO
GROUP 11. Initialization of variable or porosity fields
TKEIN=1.0;EPIN=7.5
FIINIT(U1)=UIN;FIINIT(KE)=TKEIN;FIINIT(EP)=EPIN
** Analytical solution placed in ANKE store
STORE(ANKE);REAL(AA,NN,BB,ANTKE,C2E,GX);INTEGER(JJM1)
C2E=1.92;NN=1./(1.0-C2E);AA=1./NN;BB=1.0-AA
REAL(CON1,CON2,CON3)
CON1 = TKEIN**AA ; CON2 = AA*EPIN/(UIN*TKEIN**BB)
CON3 = XULAST/2.0 ; CON2=CON2 * CON3
DO JJ=1,NX
+ PATCH(IN:JJ:,INIVAL,JJ,JJ,1,NY,1,NZ,1,1)
+ GX=XFRAC(JJ)
+ IF(JJ.NE.1) THEN
+   JJM1=JJ-1; GX=(XFRAC(JJ) + XFRAC(JJM1))
+ ENDIF
+ ANTKE=(CON1 - CON2*GX)**NN
+ INIT(IN:JJ:,ANKE,ZERO,ANTKE)
ENDDO

GROUP 13. Boundary conditions and special sources
** Inlet Boundary
PATCH(INLET,WEST,1,1,1,1,1,1,1,1)
COVAL(INLET,KE,RHO1*UIN,TKEIN);COVAL(INLET,EP,RHO1*UIN,EPIN)
GROUP 16. Termination of iterations
LSWEEP=1000;SELREF=T;RESFAC=0.0001
GROUP 21. Print-out of variables
OUTPUT(U1,N,N,N,N,N,N);OUTPUT(VIST,N,N,N,N,N,N)
GROUP 22. Monitor print-out
IXMON=NX
GROUP 23. Field print-out and plot control
NXPRIN=nx/5;ORSIZ=0.4
PATCH(PROFILE,PROFIL,1,NX,1,1,1,1,1,1)
PLOT(PROFILE,KE,0.0,1.5*TKEIN);PLOT(PROFILE,EP,0.0,2.0*EPIN)
PLOT(PROFILE,VIST,0.0,0.0)
**END OF LIBRARY CASE 230
```