Encyclopaedia Index

### 3.1 An example of time-step specification

Library case Z619

Here, for one dependent variable, H1, and four independent variables, X, Y, Z and time, PLANT is used for specification of two time-step sizes.

The first step is supposed to be equal to 5 sec, while the size of the second time step equals the smallest cell size divided by the largest velocity.

The actual solution for H1 is of no significance, so that no initial and boundary conditions for H1 are introduced. The velocities are stored, and initialised, but not solved.

The non-uniform grid and velocity field are initialised for expected size of the second time step to be 0.1 sec.

### What the user puts into Q1

```
GROUP 1. Run title and other preliminaries
TEXT(time-step CALCULATIONS
GROUP 2. Transience; time-step specification
LSTEP=2
TLAST=GRND
```

The above three lines will instruct EARTH to perform 2 time steps, the size of which being determined by the settings PLANTed in group 2.

```
** Ask PLANT to introduce first time step, 5 seconds.
{SCTS01} DT=5.
REGION(1,1,1,1,1,1,1,1)
```

The REGION qualifier is used to control this simple setting. Only the last two arguments are relevant, because the time step is necessarily the same for all cells in the domain.

```
** Find the smallest cell size at the start of each iz-slab
for the last sweep of the first time step, and place it in
the variable RG(1).
RG(1)=GREAT
DXG2D - cell size in X-direction,
DYG2D - cell size in Y-direction,
DZWNZ - cell size in Z-direction,
{SC0301} RG(1)=AMIN1(RG(1),DXG2D,DYG2D,DZWNZ)
IF(ISTEP.EQ.1.AND.ISWEEP.EQ.LSWEEP)
```

The qualifier IF restricts the whole-domain-default extents of the PLANT statement.

The arguments of IF instruct PLANT to introduce the logical conditions of the bracketed expression.

```
** Output to the globcalc of the smallest cell size using
summation in one cell at 1st time step, last sweep. This is done
for checking purposes.
{SC0302} SIZMIN=SUM(RG(1))
TEXT(Smallest cell size)
REGION(1,1,1,1,1,1,1,1) /ISWEEP.EQ.LSWEEP
```

The function SUM is used in above three lines just to dump into globcalc file the value of RG(1) headed by the character arguments of TEXT command.

```
** Find the largest velocity value just before the
second time-step calculation and place it in RG(2)
{SCTS02} RG(2)=AMAX1(RG(2),U1,V1,W1)
IF(ISTEP.EQ.2)

** Output the largest velocity using summation in one
cell at 2nd time step for checking.
{SCTS03} VELMAX=SUM(RG(2))
TEXT(Largest velocity)
REGION(1,1,1,1,1,1,2,2)
```

Auxiliary variable, RG(2), will be the largest velocity value as result of the above statements at the second time moment. It will be dumped in globcalc file and appropriately headed.

```
{SCTS04} DT=SIZMIN/VELMAX
REGION(1,1,1,1,1,1,2,2)
```

The final stage of time-step settings: second time step size, DT, is set as division of smallest size by largest velocity.

```
** The simple non-uniform grid with smallest cell size of 1 m.
GROUP 3. X-direction grid specification
GRDPWR(X,2,2.,1.)
GROUP 4. Y-direction grid specification
GRDPWR(Y,2,4.,1.)
GROUP 5. Z-direction grid specification
GRDPWR(Z,2,6.,1.)
GROUP 7. Variables stored, solved & named
SOLVE(H1)
STORE(U1,V1,W1)
GROUP 11. Initialization of variable or porosity fields
** The simple non-uniform velocity field with maximum
value, 10.0 m/s, at east-south-high corner of the domain.
FIINIT(U1)=1.;FIINIT(V1)=0.5;FIINIT(W1)=0.1
PATCH(MAXVEL,INIVAL,NX-1,NX-1,NY,NY,NZ,NZ,1,1)
COVAL(MAXVEL,U1,0.0,10.)
GROUP 15. Termination of sweeps
LSWEEP=2
GROUP 19. Data communicated by satellite to GROUND
NAMSAT=MOSG
STOP
```

### Extract from the result file

************************************************************ TIME STP= 1 SWEEP NO= 2 ZSLAB NO= 1 ITERN NO= 1 TIME = 5.000E+00 DT = 5.000E+00 ************************************************************ TIME STP= 2 SWEEP NO= 2 ZSLAB NO= 1 ITERN NO= 1 TIME = 5.100E+00 DT = 1.000E-01 ************************************************************

### What PLANT puts into GROUND

```

COMMON/SPLVRB/SIZMIN,VELMAX

CALL MAKE(DXG2D )
CALL MAKE(DYG2D )

C      Special calls name: SCTS01
IF(ISTEP.GE.1       .AND.ISTEP.LE.1       ) THEN
DO 20001 IZZ=1       ,1
DO 20001 IX=1       ,1
DO 20001 IY=1       ,1
20001 DT=5.
ENDIF
C      Special calls name: SCTS02
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISTEP.EQ.2) THEN
DO 20002 IZZ=1       ,NZ
LFU1  =ABS(ANYZ(INAME('U1    '),IZZ))
LFV1  =ABS(ANYZ(V1    ,IZZ))
LFW1  =ABS(ANYZ(W1    ,IZZ))
DO 20002 IX=1       ,NX
DO 20002 IY=1       ,NY
L0U1  =LFU1  +I
L0V1  =LFV1  +I
L0W1  =LFW1  +I
20002 RG(2)=AMAX1(RG(2),F(L0U1),F(L0V1),F(L0W1))
ENDIF
ENDIF
c     Special calls name: SCTS03
IF(ISTEP.GE.2       .AND.ISTEP.LE.2       ) THEN
VELMAX=0
DO 20003 IZZ=1       ,1
DO 20003 IX=1       ,1
DO 20003 IY=1       ,1
20003 VELMAX=VELMAX+RG(2)
ENDIF
c     Special calls name: SCTS04
IF(ISTEP.GE.2       .AND.ISTEP.LE.2       ) THEN
DO 20004 IZZ=1       ,1
DO 20004 IX=1       ,1
DO 20004 IY=1       ,1
20004 DT=SIZMIN/VELMAX
ENDIF

c     Special calls name: SC0301
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISTEP.EQ.1.AND.ISWEEP.EQ.LSWEEP) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFDXG2=L0F(DXG2D )
LFDYG2=L0F(DYG2D )
L0DZWN=L0F(DZWNZ )+IZ
DO 19301 IX=1       ,NX
DO 19301 IY=1       ,NY
L0DXG2=LFDXG2+I
L0DYG2=LFDYG2+I
19301 RG(1)=AMIN1(RG(1),F(L0DXG2),F(L0DYG2),F(L0DZWN))
ENDIF
ENDIF
ENDIF
c     Special calls name: SC0302
IF(ISTEP.GE.1       .AND.ISTEP.LE.1       ) THEN
IF(ISWEEP.EQ.LSWEEP) THEN
IF(IZ.GE.1       .AND.IZ.LE.1       ) THEN
IF(IZ.EQ.1       ) SIZMIN=0
DO 19302 IX=1       ,1
DO 19302 IY=1       ,1
19302 SIZMIN=SIZMIN+RG(1)
ENDIF
ENDIF
ENDIF

LU=1
OPEN(UNIT=LU,FILE='GLOBCALC',STATUS='UNKNOWN')
WRITE(LU,*) 'Global calculations:'
WRITE(LU,*) 'SMALLEST CELL SIZE            ',' = ',SIZMIN
WRITE(LU,*) 'LARGEST VELOCITY              ',' = ',VELMAX
CLOSE(LU)

```

### 3.2 An example of expanding and contracting grid

Library case Z622

In this example, attention is focussed wholly on how to use PLANT for expanding and/or contracting the Y-extent of the grid.

It is useful in parabolic calculations when PARAB=T. The grids created this way can also be used for independent purposes.

The example is made for Y-direction grid specification. The PLANT statements headed by SCXS?? should perform an exactly corresponding functions for the X-direction domain width, XULAST, to that performed here for YVLAST, provided the pointer AZXU is set to GRND.

### What the user puts into Q1

To generate the above grid user should put the following in the Q1 file:

```
AZYV=GRND
NAMSAT=MOSG

REGION(1,1,1,1) /IZ.LE.20

REGION(1,1,1,1) /IZ.GT.20.AND.IZ.LE.40

REGION(1,1,1,1) /IZ.GT.40
```

The above three statements make the domain width a power-law function of the downstream distance of the current slab (ZW).

The power is changed with Z-direction in step-wise manner. It is positive over first 20 slabs, negative over next 20 slabs and positive again for the remainder of the domain.

If the z-direction step size is made proportional to YVLAST by setting AZDZ=PROPY and DZW1 as follows

NZ=100; AZDZ=PROPY; DZW1=0.05

the resulting grid will be as shown on the picture.

### What PLANT puts into GROUND

In respond of settings in Q1 the following Ground coding is created by PLANT:-

```
c     Special calls name: SCYS01
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.LE.20) THEN
IF(IZ.GE.1       ) THEN
DO 40001 IX=1       ,1
DO 40001 IY=1       ,1
ENDIF
ENDIF
ENDIF
c     Special calls name: SCYS02
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GT.20.AND.IZ.LE.40) THEN
IF(IZ.GE.1       ) THEN
DO 40002 IX=1       ,1
DO 40002 IY=1       ,1
ENDIF
ENDIF
ENDIF
c     Special calls name: SCYS03
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GT.40) THEN
IF(IZ.GE.1       ) THEN
DO 40003 IX=1       ,1
DO 40003 IY=1       ,1
ENDIF
ENDIF
ENDIF
```

### 3.3 An example of adaptive Z-direction step size

Library case Z623

The case is concerned with the calculation of z-direction step size for parabolic problems. The case 150 is loaded from the core library and appended by PLANT blocks.

They create the GROUND codings for composite function governing the distribution of grid steps in Z-direction: DZ is proportional to the Y-extent of the grid for the first two steps and each next step size is reduced as longitudinal centre-line temperature decay dictates.

The consequence is that the Z-direction step size gradually reduces towards the domain end as shown.

The example also demonstrates the ability to overwrite some PIL fragments by the PLANT sequences.

### What the user puts into Q1

The user needs to put the following in the Q1 file

```
GROUP 5. Z-grid
AZDZ=GRND
{SCZS01} DZ=0.1*YVLAST
IF(IZ.LE.2)
{SCZS02} DZ=DZL*(1.-0.02*TEMP[,1,-1]/:TJET:)
IF(IZ.GT.2)
```

The above two PLANT blocks overwrite the correspondong actions in PIL of Library Case 150.

They make the step size, DZ, a fixed fraction of the calculation domain width, YVLAST, for the steps less or equal 2. If the step number is greater than 2, the each next step size is made a fraction of the thickness for previous z-slab, DZL.

The fraction factor is supposed to be adaptive to the upstream (minus unity third argument in square bracets) temperature, TEMP, at the first cell in Y-direction ( unity second argument in square brackets).

The calculations result in the following grid.

### What PLANT puts into GROUND

```
c     Special calls name: SCZS01
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.LE.2) THEN
IF(IZ.GE.1       ) THEN
DO 50001 IX=1       ,NX
DO 50001 IY=1       ,NY
50001 DZ=0.1*YVLAST
CALL FN1(LXYDZ,DZ)
ENDIF
ENDIF
ENDIF
c     Special calls name: SCZS02
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GT.2) THEN
IF(IZ.GE.1       ) THEN
LFTEMP=L0F(INAME('TEMP  '))
DO 50002 IX=1       ,NX
DO 50002 IY=1       ,NY
L0TEMP=LFTEMP+I
50002 DZ=DZL*(1.-0.02*F(L0TEMP+1-IY-NFM)/1)
CALL FN1(LXYDZ,DZ)
ENDIF
ENDIF
ENDIF
```

### 3.4 An example of analytical BFC grid generation

Library case Z607 Library case Z608 Library case Z609 Library case Z501 Library case Z502 Library case Z503 Library case Z504 Library case Z505 Library case Z506 Library case Z507 Library case Z508 Library case Z509 Library case Z510 Library case Z511 Library case Z512 Library case Z513 Library case Z514

PLANT can be used as a formula-based BFC grid generator. The complex BFC grids are created in response of typing in Q1 parameterised analytical expressions as exemplified below.

### What the user puts into Q1

```
REAL(LENGTH,TWOPI,LITTLER,PARAM)
LENGTH=10.0; LITTLER=1.0;TWOPI=2.0*3.14157
PARAM=1.0
{MXYZ01}  XC=ABS(:PARAM:*COS(:LENGTH:*FLOAT(K-1)/FLOAT(NZ)))+\$
:LITTLER:*FLOAT(J-1)/FLOAT(NY)*\$
COS(:TWOPI:*FLOAT(I-1)/FLOAT(NX))
{MXYZ01}  YC=ABS(:PARAM:*COS(:LENGTH:*FLOAT(K-1)/FLOAT(NZ)))+\$
:LITTLER:*FLOAT(J-1)/FLOAT(NY)*\$
SIN(:TWOPI:*FLOAT(I-1)/FLOAT(NX))
{MXYZ01}  ZC=:LENGTH:*FLOAT(K-1)/FLOAT(NZ)
IF(ISTEP.EQ.1.AND.ISWEEP.EQ.1)
```
The above three statements followed by IF command provide the calculation of cartesian coordinates for cell corners of the grid fitted the corrugated, in accord with abs(cosZ) function, circular pipe of 1m radius and 10m length.

The grid is uniform in both direction. The generation is made at the first sweep of the first time step.

The calculation results in the grid shown on the next picture.

The formula parameters in above setttings are the pipe length, LENGTH, its radius, LITTLER, and corrugation factor, PARAM.

If some of them are set to be, say, a function of time the moving BFC grid can be esily generated.

For example, the replacement of :PARAM: by :PARAM:*TIM would result in the number of BFC grid generated for each time moment.

The next three pictures illustrate the grids for suuccesive time moments as follows:

Time = 0.0 sec

Time = 0.25 sec

Time = 0.75 sec

### What PLANT puts into GROUND

```

c    Special calls name: MXYZ01
IF(BFC) THEN
IF(ISTEP.EQ.1.AND.ISWEEP.EQ.1)THEN
DO 19201 K=1       ,NZ      +1
DO 19201 I=1       ,NX      +1
DO 19201 J=1       ,NY      +1
XC=ABS(7.5000E-01*COS(10*FLOAT(K-1)/FLOAT(NZ)))+1*FLOAT
1(J-1)/FLOAT(NY)*COS(6.2831E+00*FLOAT(I-1)/FLOAT(NX))
YC=ABS(7.5000E-01*COS(10*FLOAT(K-1)/FLOAT(NZ)))+1*FLOAT
1(J-1)/FLOAT(NY)*SIN(6.2831E+00*FLOAT(I-1)/FLOAT(NX))
ZC=10*FLOAT(K-1)/FLOAT(NZ)
CALL SECRNS(XCORNR,YCORNR,ZCORNR,I,J,K,XC,YC,ZC)
19201  CONTINUE
CALL BGEOM(1)
CALL BGEOM(2)
CALL DUMPS(CSG1(1:1),CSG2(1:1),0,1,0,0)
ENDIF
ENDIF
```

### 3.5 An example of convection-flux alteration

Library case Z614 Library case Z621

This calculation is of scalar source dispersion in 45 degree uniform flow. PLANt is used to alter the convection fluxes so that the plume can propagate in the direction opposite to the flow.

The result vector field and the dispersion plume are as shown.

The case presented exploits the ability of PLANT for easy indicial expression operations to create multi-domain computational space.

### What the user puts into Q1

Although the problem in question is 2D, NX*NY=20*20, in nature the provision is made below for cartesian box to have 2 slabs.

```
GROUP 3. X-direction grid specification
GRDPWR(x,20,20.,1.0)
GROUP 4. Y-direction grid specification
GRDPWR(Y,20,20.,1.0)
GROUP 5. Z-direction grid specification
GRDPWR(Z,2,2.,1.0)
```
The convection-diffusion transport of scalar in prescribed velocity field will be considered here, so that:

```
GROUP 7. Variables stored, solved & named
SOLVE(H1)
STORE(U1,V1,W1,HPOR)
GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,Y,Y,Y,Y,Y)
```
The nullification, by two commands below, of high face porosities provides the independency of the slab sub-domains.

```
GROUP 11. Initialization of variable or porosity fields
FIINIT(HPOR)=0.0
```

The following set of initialisations make the 45 degree flow of 5 m/s from south-east edge of the domain. It will be maintained as 1st and 2nd slab velocity fields.

FIINIT(V1)=5.0;FIINIT(U1)=-5.0

Although, the velocity field at the second slab is the same as for first one, the add-extra-velocity option is activated as PLANT settings below tell.

```

NAMSAT=MOSG
REGION(,NX-1) /IZ.EQ.2
REGION(,,,NY-1)
IF(IZ.EQ.2)
```

The extra velocities added to the flow velocity components alters the convection fluxes to become opposite to ones at first slab.

Please note the differences in the REGION qualifier. They are attributed to the staggered nature of velocity nodes. The alternative use of either switch, /IZ.EQ.2, or IF command limits the Z-direction extent of velocity alterations.

The expected distribution of convected property H1 fixed at the middle of the second slab by

```
PATCH(FIXSOR,CELL,nx/2,nx/2,NY/2,NY/2,2,2,1,1)
COVAL(FIXSOR,H1,FIXVAL,1.0)
```

results in the plume towards the north-west corner of the domain.

### What PLANT puts into GROUND

In the GROUND file, PLANT inserts the following FORTRAN coding,

c Special calls name: SCUF01

```
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.EQ.2) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFU1  =L0F(U1    )
DO 81001 IX=1       ,NX-1
DO 81001 IY=1       ,NY
L0VELA=LFVELA+I
L0U1  =LFU1  +I
81001 F(L0VELA  )=-2.*F(L0U1+NFM*(1-IZ))
ENDIF
ENDIF
ENDIF
c     Special calls name: SCVF01
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.EQ.2) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFV1  =L0F(V1    )
DO 83001 IX=1       ,NX
DO 83001 IY=1       ,NY-1
L0VELA=LFVELA+I
L0V1  =LFV1  +I
83001 F(L0VELA  )=-2.*F(L0V1+NFM*(1-IZ))
ENDIF
ENDIF
ENDIF
```

### 3.6 An example of a non-linear property correlation

Library case Z612 The case presented is entitled 'Lid-driven flow with property variations'.

Briefly, this calculation is of flow in a cavity, the top wall of which moves with constant velocity. The top wall is at one temperature, the moving wall is at a different temperature and the other walls are adiabatic.

The calculation is made for Reynolds number Re = 100 for Prandtl number Pr = 1 and for an exponent type viscosity formula of the type,

EMU = EMU0 * EXP(-beta.(T - T0))

The geometry and vector and temperature-contour results are as shown.

### What the user puts into Q1

```
GROUP 9. Properties of the medium (or media)
PRNDTL(TEMP)=0.7;ENUL=GRND
REAL(ENULR,EXPO,TEMPR)
ENULR=1.E-3
EXPO=-1.6
TEMPR=0.0
{PRPT01} VISL=:ENULR:*EXP(:EXPO:*(TEMP-:TEMPR:))
```

The result viscosity field is highly non-uniform.

### What PLANT puts into GROUND

```

In the GROUND file, PLANT inserts the following FORTRAN coding,

c     Property name: PRPT01
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFTEMP=L0F(INAME('TEMP  '))
DO 96001 IX=IXF,IXL
DO 96001 IY=IYF,IYL
L0TEMP=LFTEMP+I
L0VISL=LFVISL+I
96001 F(L0VISL  )=1.0000E-03*EXP(-1.6000E+00*(F(L0TEMP)-0))
ENDIF
ENDIF
```

### 3.7 An example of interphase transport

Library case Z606

TEXT(VEHICULAR EXHAUST DISPERSION IN RAINFALL

This case simulates interactions between air and rainfall in the street between a tall and a low-rise neighbouring building.

The pollution is caused by inter-phase transport from gaseous to liquid phase.

The contaminant concentration in liquid phase is depicted on the picture.

PLANT instructions provided in Group 10 make the special GROUND codings.

### What the user puts into Q1

```
*  GROUP 10. Inter-phase-transfer processes and properties.
** Set a constant inter-phase friction coefficient.

CFIPS=gravac*(rho2-rho1)/fallvl
** Instruct PLANT to code the inter-phase  transfer  coefficients
and PHI differences between the phases.

CINT(H1)=GRND
{PRPT01} COI1(H1)=10.*MASS2                         (Equation 1)
CINT(H2)=1.e10
PHINT(H1)=GRND
{PRPT02} FII1(H1)=H2                                (Equation 2)
PHINT(H2)=GRND
{PRPT03} FII2(H2)=H1                                (Equation 3)

```

### What PLANT puts into GROUND

```

For Equation 2, the following Ground coding was created by PLANT:-

c     Property name: PRPT02
IF(INDVAR.EQ.INAME('H1    ')) THEN
LFFII1 =L0F(FII1  )
LFH2  =L0F(H2    )
DO 98002 IX=IXF,IXL
DO 98002 IY=IYF,IYL
L0H2  =LFH2  +I
L0FII1=LFFII1+I
98002 F(L0FII1  )=F(L0H2)
ENDIF

For Equation 3, the following Ground coding was created by PLANT:-

c     Property name: PRPT03

IF(INDVAR.EQ.INAME('H2    ')) THEN
LFFII2 =L0F(FII2  )
LFH1  =L0F(H1    )
DO 99003 IX=IXF,IXL
DO 99003 IY=IYF,IYL
L0H1  =LFH1  +I
L0FII2=LFFII2+I
99003 F(L0FII2  )=F(L0H1)
ENDIF

For Equation 1, the following Ground coding was created by PLANT:-

c  Property name: PRPT01

IF(INDVAR.EQ.INAME('H1    ')) THEN
LFCOI1 =L0F(COI1  )
LFMAS2=L0F(MASS2 )
DO 10301 IX=IXF,IXL
DO 10301 IY=IYF,IYL
L0MAS2=LFMAS2+I
L0COI1=LFCOI1+I
10301 F(L0COI1  )=10.*F(L0MAS2)
ENDIF
```

The reader is asked to compare these latter three panels with what was placed in the Q1 file.

### 3.8 An example of devices for setting initial and porosity fields

Library case Z605

The settings which appear below illustrate some of the existing initialisation techniques available in PLANT.

They cover:

### What the user puts into Q1

a) For flow field initialistions by parametric analytics

```
PATCH(INI1,INIVAL,1,NX,1,NY,1,1,1,1)
{INIT01} VAL=XU2D-10.
COVAL(INI1,U1,zero,GRND)
{INIT02} VAL=-(YV2D-10.)
COVAL(INI1,V1,zero,GRND)
```

By above settings the first slab of 3D domain is initialized by stagnation point flow with the cartesian components as follows: U1 = X - 10 and V1 = 10- Y.

Resulting initial flow field is as shown.

The settings below initialize the flow field in 2nd slab by solid body rotation components, as follows:

```
PATCH(INI2,INIVAL,1,NX,1,NY,2,2,1,1)
{INIT03} VAL=YG2D-10.
COVAL(INI2,U1,zero,GRND)
{INIT04} VAL=-(XG2D-10.)
COVAL(INI2,V1,zero,GRND)
```

Resulting initialization is as shown.

b) For manipulating with initial fields

```
PATCH(INI3,INIVAL,1,NX,1,NY,3,3,1,1)
{INIT05} VAL=-U1[,,-1]+U1[,,1]
COVAL(INI3,U1,zero,GRND)
{INIT06} VAL=-V1[,,2]+V1[,,-2]
COVAL(INI3,V1,zero,GRND)
```

By above settings, the 3rd slab of 3D domain is initialized by superposition of velocities at 2nd and 1st slabs.

The composite initialization is as shown.

c) For MARKing sub-domains to create arbitrary initial fields

```
FIINIT(MARK)=1.0
PATCH(INIT70,INIVAL,1,NX/2,1,NY,4,4,1,1)
{INIT70} VAL=XYELLP(2.,10.,10.,8.,8.,0.,0.)
INIT (INIT70,MARK,0.,GRND)
```

In above statement, XYELLP function is used to make the half-circle of 16 m diameter as follows:

• i) In the west half of 4th slab,
• ii) place the ellipse marked by 2 (1st argument), with the centre at XC=10 m (2nd argument) and YC=10 m (3rd argument) and both half-axes equal to 8 m ( 4th and 5th arguments).
• iii) The 6th and 7th function arguments are insignificant for the circle shape and CARTES=T.

```
{SC0301} U1=U1[,,2]
REGION() 2 /ISWEEP.EQ.1

{SC0302} V1=V1[,,2]
REGION() 2 /ISWEEP.EQ.1
```

The above two statements followed by the REGION qualifiers with parameters and switches set the velocity components of 4th slab region marked by MARK=2 to be equal of those of 2nd slab. It is done at the start of IZ=2 slab for the first sweep only.

By this way arbitrary complex flow field can be initialized as picture shows.

The same technique provides the possibility to distribute the PRPS values to describe the complex geometry on Cartesian, polar and BFC grids by blocking-off.

Thefirst picture shows how the resulting shape may look like on Cartesian grid.

The next picture depicts the blocking-off polar grid.

Finally, the result of blocking-off the BFC grid is shown.

PLANT is equipped with a number of geometrical functions which can be used in Q1 settings to make all above actions easy and straightforward.

### 3.9 An example of non-linear boundary conditions and source laws

Library case Z6106

TEXT(3D STEADY HEAT CONDUCTION IN A CUBE RADIATION AND CONVECTION FROM A HOT BLOCK

This example exemplifies the use of PLANT for external heat-transfer laws.

This feature allows inclusion of external heat loss by way of radiation using the distribution of radiosity, i.e., (sigma(Text**4-Ts**4)) and forced convection a(Text-Ts).

The geometry and temperature distribution are shown on the picture.

### What the user puts into Q1

```

GROUP 13. Boundary conditions and special sources

{SORC01} CO=RG(1)*(RG(2)**2+H1**2)*(RG(2)+H1)
{SORC01} VAL=RG(2)
```

### What PLANT puts into GROUND

```
c     Source name: SORC01
LFCO  =L0F(CO )
LFH1  =L0F(H1    )
DO 13701 IX=IXF,IXL
DO 13701 IY=IYF,IYL
L0H1  =LFH1  +I
13701 F(LFCO +I)=RG(1)*(RG(2)**2+F(L0H1)**2)*(RG(2)+F(L0H1))
ENDIF
LFVAL =L0F(VAL)
DO 13801 IX=IXF,IXL
DO 13801 IY=IYF,IYL
13801 F(LFVAL+I)=RG(2)
ENDIF
```

### 3.10 An example of setting residual reference values

Click to return here after viewing pictures

The example is concerned with the calculation of reference residuals for dependent variable G as an arithmetic mean of its net generation rate over the slab.

### What the user puts into Q1

```

{SC0655} RES=SUM(VOL*(GENG-2.*:RHO1:*EPKE*G)/(NY*NZ))
{SC0656} RESREF(G)=RES
```

By above two statements the reference residuals for G is calculated at the end of each z-slab.

### What PLANT puts into GROUND

```

COMMON/SPLVRB/RES

c     Special calls name: SC0655
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
IF(IZ.EQ.1       ) RES   =0
LFVOL =L0F(VOL   )
LFGENG=L0F(INAME('GENG  '))
LFEPKE=L0F(INAME('EPKE  '))
LFG   =L0F(INAME('G     '))
DO 19655 IX=1       ,NX
DO 19655 IY=1       ,NY
L0VOL =LFVOL +I
L0GENG=LFGENG+I
L0EPKE=LFEPKE+I
L0G   =LFG   +I
19655 RES=RES+F(L0VOL)*(F(L0GENG)-2.*1*F(L0EPKE)*F(L0G))/
1(NY*NZ)
ENDIF
ENDIF
c     Special calls name: SC0656
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
DO 19656 IX=1       ,NX
DO 19656 IY=1       ,NY
19656 RESREF(INAME('G'))=RES
ENDIF
ENDIF

c  * ------------------- SECTION 8 ---- Finish of time step.
LU=1
OPEN(UNIT=LU,FILE='GLOBCALC',STATUS='UNKNOWN')
WRITE(LU,*) 'Global calculations:'
WRITE(LU,*) 'RES                           ',' = ',RES
CLOSE(LU)
```

### 3.11 An example of setting under-relaxation devices

Library case Z613 Click to return here after viewing pictures

Two types of self-steering false-time under-relaxations for the velocities are introduced by way of PLANT, namely global and local self-steering.

To see the power of relaxation , the solution process is divided into 3 stages:

1. No relaxation for ISWEEP less and equal 100,
2. Global relaxation for ISWEEP larger than 100 and less or equal 200 and
3. Local self-steering relaxation for ISWEEP larger than 200.

The picture shows the convergence process which is clearly divided into three stages described above.

```
{SC0201} RG(2)=AMIN1(XULAST/FLOAT(NX),YVLAST/FLOAT(NY),\$
ZWLAST/FLOAT(NZ))/\$
AMAX1(U1,:FLO1:/2.)
REGION(1,1,2,3,2,2)
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200)

{SC0202} DTFALS(U1)=RG(2)
REGION(1,1,1,1,1,1)
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200)
```

In the above, global under-relaxation is introduced by PLANTed codings for DTFALS(U1) at the start of each sweep.

It is assumed to be equal to the smallest of the cell sizes divided by the largest of inlet mass flux velocity and local velocity magnitude normal to the inlet plane.

It is applied over the whole domain for the velocity in question IF isweep is greater than 100 but less or equal than 200.

Here and for next two statemnts, command REGION with unity arguments is used to economize the operations needed for equivalences.

```
{SC0203} DTFALS(V1)=RG(2)
REGION(1,1,1,1,1,1)
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200)
```

The above settings do for DTFALS(V1) what has been done for DTFALS(U1) above.

```
{SC0204} DTFALS(W1)=RG(2)
REGION(1,1,1,1,1,1)
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200)
```

The above settings do for DTFALS(W1) what has been done for DTFALS(U1).

```
PATCH(RELAX,PHASEM,1,NX,1,NY,1,NZ,1,1)
{SORC97} CO=1./TFAL
COVAL(RELAX,U1,GRND,SAME)
IF(ISWEEP.GT.200)
{SORC98} CO=1./TFAL
COVAL(RELAX,V1,GRND,SAME)
IF(ISWEEP.GT.200)
{SORC99} CO=1./TFAL
COVAL(RELAX,W1,GRND,SAME)
IF(ISWEEP.GT.200)
```

In the above settings, local self-steering under-relaxation is introduced through the sources of momentum for the whole domain defined by PATCH named RELAX, for which:

• TYPE is PHASEM,
• VALue is SAME,
• COefficient, is set to reciprocal of false-time step.

It is applied for each sweep greater than 200.

```
STORE(TFAL);OUTPUT(TFAL,Y,Y,Y,Y,Y,Y)
{SC0301} TFAL=1/(SQRT(U1**2+W1**2+V1**2)/\$
AMIN1(DXU2D*1,AMIN1(DYV2D*1,DZ*1))+\$
RG(1)/AMIN1(DXU2D*1,AMIN1(DYV2D*1,DZ*1))**2)
IF(ISWEEP.GT.200)
```

The reciprocal of local self-steering false-time step is set to the local velocity vector magnitude divided by smallest distance between walls of continuity cell in question plus local diffusivities, i.e. kinematic viscosities, divided by the smallest distance squarred.

The variable TFAL, false-time, is provided to assist the computations. It is calculated right at the start of each IZ-slab for all sweeps greter than 200 and can be used to monitor the variation of local magnitudes of false-time steps.

### What PLANT puts into GROUND

```
C*The following codings were created by PLANT in the GROUND file:
c     Source name: SORC97
IF(INDVAR.EQ.INAME('U1    ').AND.NPATCH.EQ.'RELAX   ') THEN
IF(ISWEEP.GT.200) THEN
LFCO  =L0F(CO )
LFTFAL=L0F(INAME('TFAL  '))
DO 13797 IX=IXF     ,IXL
DO 13797 IY=IYF     ,IYL
L0TFAL=LFTFAL+I
13797 F(LFCO +I)=1./F(L0TFAL)
ELSE
LFCO =L0F(CO)
DO IX=IXF,IXL
DO IY=IYF,IYL
F(LFCO +I)=0.
ENDDO
ENDDO
ENDIF
ENDIF
c     Source name: SORC98
IF(INDVAR.EQ.INAME('V1    ').AND.NPATCH.EQ.'RELAX   ') THEN
IF(ISWEEP.GT.200) THEN
LFCO  =L0F(CO )
LFTFAL=L0F(INAME('TFAL  '))
DO 13798 IX=IXF     ,IXL
DO 13798 IY=IYF     ,IYL
L0TFAL=LFTFAL+I
13798 F(LFCO +I)=1./F(L0TFAL)
ELSE
LFCO =L0F(CO)
DO IX=IXF,IXL
DO IY=IYF,IYL
F(LFCO +I)=0.
ENDDO
ENDDO
ENDIF
ENDIF
c     Source name: SORC99
IF(INDVAR.EQ.INAME('W1    ').AND.NPATCH.EQ.'RELAX   ') THEN
IF(ISWEEP.GT.200) THEN
LFCO  =L0F(CO )
LFTFAL=L0F(INAME('TFAL  '))
DO 13799 IX=IXF     ,IXL
DO 13799 IY=IYF     ,IYL
L0TFAL=LFTFAL+I
13799 F(LFCO +I)=1./F(L0TFAL)
ELSE
LFCO =L0F(CO)
DO IX=IXF,IXL
DO IY=IYF,IYL
F(LFCO +I)=0.
ENDDO
ENDDO
ENDIF
ENDIF
CALL MAKE(DXU2D )
CALL MAKE(DYV2D )
c     Special calls name: SC0201
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200) THEN
DO 19201 IZZ=2       ,2
LFU1  =ABS(ANYZ(U1    ,IZZ))
DO 19201 IX=1       ,1
DO 19201 IY=2       ,3
L0U1  =LFU1  +I
19201 RG(2)=AMIN1(XULAST/FLOAT(NX),YVLAST/FLOAT(NY),ZWLAST/FL
1OAT(NZ))/AMAX1(F(L0U1),1.0000E-01/2.)
ENDIF
ENDIF
c     Special calls name: SC0202
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200) THEN
DO 19202 IZZ=1       ,1
DO 19202 IX=1       ,1
DO 19202 IY=1       ,1
19202 DTFALS(U1)=RG(2)
ENDIF
ENDIF
c     Special calls name: SC0203
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200) THEN
DO 19203 IZZ=1       ,1
DO 19203 IX=1       ,1
DO 19203 IY=1       ,1
19203 DTFALS(V1)=RG(2)
ENDIF
ENDIF
c     Special calls name: SC0204
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISWEEP.GT.100.AND.ISWEEP.LE.200) THEN
DO 19204 IZZ=1       ,1
DO 19204 IX=1       ,1
DO 19204 IY=1       ,1
19204 DTFALS(W1)=RG(2)
ENDIF
ENDIF
c     Special calls name: SC0301
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(ISWEEP.GT.200) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFTFAL =L0F(INAME('TFAL  '))
LFU1  =L0F(U1    )
LFW1  =L0F(W1    )
LFV1  =L0F(V1    )
LFDXU2=L0F(DXU2D )
LFDYV2=L0F(DYV2D )
DO 19301 IX=1       ,NX
DO 19301 IY=1       ,NY
L0TFAL=LFTFAL+I
L0U1  =LFU1  +I
L0W1  =LFW1  +I
L0V1  =LFV1  +I
L0DXU2=LFDXU2+I
L0DYV2=LFDYV2+I
19301 F(L0TFAL  )=1/(SQRT(F(L0U1)**2+F(L0W1)**2+F(L0V1)*
1*2)/AMIN1(F(L0DXU2)*1,AMIN1(F(L0DYV2)*1,DZ*1))+RG(1)/AM
1IN1(F(L0DXU2)*1,AMIN1(F(L0DYV2)*1,DZ*1))**2)
ENDIF
ENDIF
ENDIF
```

### 3.12 An example of setting limits and increments

Click to return here after viewing pictures Library case Z620

The limit on variable value, TEMP, to be the largest stored/solved field variable value, TAD, is introduced here by way of two-stage settings.

### What the user puts into Q1

```

{SC0302} VARMAX(TEMP)=RG(1)
```

The two statements above provide that the value of VARMAX(TEMP) is set equal to RG(1) which is the largest of the TAD values.

It is done at the start of each iz-slab.

### What PLANT puts into GROUND

```

c     Special calls name: SC0301
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
DO 19302 IX=1       ,NX
DO 19302 IY=1       ,NY
ENDIF
ENDIF
c     Special calls name: SC0302
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
DO 19303 IX=1       ,NX
DO 19303 IY=1       ,NY
19302 VARMAX(INAME('TEMP'))=RG(1)
ENDIF
ENDIF
```

### 3.13 An example of output control

Library case Z611 Library case Z615 Library case Z616 Library case Z350 Library case Z450 Library case Z451 Library case Z452 Library case Z453 Library case Z454 Library case Z455 Library case Z456 Library case Z457 Library case Z458 Library case Z459 Library case Z147

The case presented is concerned with the calculation of the temperature distribution in the packed bed of melting solids.

The output problem is to re-calculate the final temperature field to get an estimation of mass fraction which remains in solid state.

The solid fraction is supposed to be determined from the following equation:

FS = ( Liquidus temp - Temp ) / (Liquidus temp - Solidus temp )

The bottom picture of the panel shows the distribution of temperature ( red - 2478 K, blue - 378 K ).

The upper picture represents the field of solid mass fraction.

To get required output user needs to put the following in the Q1 file

```
REAL(TLIQ,TSOL)
TLIQ=1750.
TSOL=1310.
STORE(FS);FIINIT(FS)=0.0
{SC0601} FS=AMIN1(1.,AMAX1((:TLIQ:-TEMP)/(:TLIQ:-:TSOL:),0.0))
```

### What PLANT puts into GROUND

```

c     Special calls name: SC0601
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFFS   =L0F(INAME('FS    '))
LFTEMP  =L0F(INAME('TEMP  '))
DO 19601 IX=1       ,NX
DO 19601 IY=1       ,NY
L0FS  =LFFS  +I
L0TEMP  =LFTEMP  +I
F(L0FS    )=AMIN1(1.,AMAX1((1750-F(L0TEMP))/(1750-1310
1),0.0))
19601  CONTINUE
ENDIF
ENDIF
```

### 3.14 Confined turbulent flame

Click to return here after viewing pictures

The problem of confined flame, settings of which will be presented below, requires:

• the temperature, TEMP, and
• the mass fraction of mixture components, FUEL, OXID and PROD,
to be re-calculated in accord with then solved-for variable, MIXF, field.

The flow density should then be brought in consistency with the mixture composition, temperature and pressure.

The SCRS (ie Simple Chemically Reacting System) model of gaseous combustion provides the formulae required.

The calculations result in the following:-

* for temperature;

* for fuel mass fraction and

* for density distribution.

### What the user puts into Q1

```

RHO1=GRND
{PRPT01}DEN1=FUEL/:WFU:+OXID/:WAIR:+PROD/:WPR:
{PRPT02}DEN1= (PRESS0+P1)/(8314.*DEN1*TEMP+TINY)

{SC0661} FUEL=AMAX1(0.,:RIMF:*FMIX-:RIMF:*:FSTOI:)
{SC0662} OXID=AMAX1(0.,1.-FMIX/:FSTOI:)
{SC0664} TEMP=(H1-:HFU:*FUEL)/:CP:
{SC0665} PROD=1.-OXID-FUEL
```

### What PLANT, as a consequence, puts into GROUND

```

c    Property name: PRPT01
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFDEN1 =L0F(AUX(DEN1  ))
LFFUEL=L0F(INAME('FUEL  '))
LFOXID=L0F(INAME('OXID  '))
LFPROD=L0F(INAME('PROD  '))
DO 90101 IX=IXF     ,IXL
DO 90101 IY=IYF     ,IYL
L0DEN1=LFDEN1+I
L0FUEL=LFFUEL+I
L0OXID=LFOXID+I
L0PROD=LFPROD+I
90101 F(L0DEN1  )=F(L0FUEL)/28+F(L0OXID)/29+F(L0PROD)/38
ENDIF
ENDIF
c     Property name: PRPT02
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFDEN1 =L0F(AUX(DEN1  ))
LFP1  =L0F(P1    )
LFTEMP=L0F(INAME('TEMP  '))
DO 90102 IX=IXF     ,IXL
DO 90102 IY=IYF     ,IYL
L0DEN1=LFDEN1+I
L0P1  =LFP1  +I
L0TEMP=LFTEMP+I
90102 F(L0DEN1  )=(PRESS0+F(L0P1))/(8314.*F(L0DEN1)*F(L0TEMP)+
1TINY)
ENDIF
ENDIF
c     Special calls name: SC0601
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFFUEL =L0F(INAME('FUEL  '))
LFFMIX=L0F(INAME('FMIX  '))
DO 19601 IX=1       ,NX
DO 19601 IY=1       ,NY
L0FUEL=LFFUEL+I
L0FMIX=LFFMIX+I
19601 F(L0FUEL  )=AMAX1(0.,1.4484E+00*F(L0FMIX)-1.4484E+00*3.
10958E-01)
ENDIF
ENDIF
c     Special calls name: SC0602
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFOXID =L0F(INAME('OXID  '))
LFFMIX=L0F(INAME('FMIX  '))
DO 19602 IX=1       ,NX
DO 19602 IY=1       ,NY
L0OXID=LFOXID+I
L0FMIX=LFFMIX+I
19602 F(L0OXID  )=AMAX1(0.,1.-F(L0FMIX)/3.0958E-01)
ENDIF
ENDIF
c     Special calls name: SC0603
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFTEMP =L0F(INAME('TEMP  '))
LFH1  =L0F(H1    )
LFFUEL=L0F(INAME('FUEL  '))
DO 19603 IX=1       ,NX
DO 19603 IY=1       ,NY
L0TEMP=LFTEMP+I
L0H1  =LFH1  +I
L0FUEL=LFFUEL+I
19603 F(L0TEMP  )=(F(L0H1)-1.0000E+07*F(L0FUEL))/1500.
ENDIF
ENDIF
c     Special calls name: SC0604
IF(ISTEP.GE.1       .AND.ISTEP.LE.LSTEP   ) THEN
IF(IZ.GE.1       .AND.IZ.LE.NZ      ) THEN
LFPROD =L0F(INAME('PROD  '))
LFOXID=L0F(INAME('OXID  '))
LFFUEL=L0F(INAME('FUEL  '))
DO 19604 IX=1       ,NX