- time-step specification
- expanding and contracting grid
- adaptive grid-size calculation
- analytical BFC grid generation
- convection-fluxes alteration
- a non-linear property correlation
- interphase transport
- devices for setting initial variables and porosity fields
- non-linear boundary conditions and source laws
- setting residual-reference values
- setting under-relaxation devices
- setting limits and increments
- output control
- Confined turbulent flame
- A shell-and-tube heat exchanger

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.

GROUP 1. Run title and other preliminaries TEXT(time-step CALCULATIONS GROUP 2. Transience; time-step specification STEADY=F 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. INIADD=F 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

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 IADD=NY*(IX-1) DO 20001 IY=1 ,1 I=IY+IADD 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 IADD=NY*(IX-1) DO 20002 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 20003 IY=1 ,1 I=IY+IADD 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 IADD=NY*(IX-1) DO 20004 IY=1 ,1 I=IY+IADD 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 IADD=NY*(IX-1) DO 19301 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 19302 IY=1 ,1 I=IY+IADD 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)

Click here to return to "Contents"

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.

This is what the generated grid will look like.

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

AZYV=GRND NAMSAT=MOSG{SCYS01} YRAT=1.+:POWER:*DZ/(ZWADD+ZW) REGION(1,1,1,1) /IZ.LE.20

{SCYS02} YRAT=1.+:POWER:*(-1.)*DZ/(ZWADD+ZW) REGION(1,1,1,1) /IZ.GT.20.AND.IZ.LE.40

{SCYS03} YRAT=1.+:POWER:*1.*DZ/(ZWADD+ZW) 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.

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 IADD=NY*(IX-1) DO 40001 IY=1 ,1 I=IY+IADD 40001 YRAT=1.+5.0000E-01*DZ/(ZWADD+ZW) 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 IADD=NY*(IX-1) DO 40002 IY=1 ,1 I=IY+IADD 40002 YRAT=1.+5.0000E-01*(-1.)*DZ/(ZWADD+ZW) 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 IADD=NY*(IX-1) DO 40003 IY=1 ,1 I=IY+IADD 40003 YRAT=1.+5.0000E-01*1.*DZ/(ZWADD+ZW) ENDIF ENDIF ENDIF

Click here to return to "Contents"

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.

The user needs to put the following in the Q1 file

LOAD(150) 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.

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 IADD=NY*(IX-1) DO 50001 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 50002 IY=1 ,NY I=IY+IADD L0TEMP=LFTEMP+I 50002 DZ=DZL*(1.-0.02*F(L0TEMP+1-IY-NFM)/1) CALL FN1(LXYDZ,DZ) ENDIF ENDIF ENDIF

Click here to return to "Contents"

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.

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.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 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

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

Click here to return to "Contents"

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.

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.

The convection-diffusion transport of scalar in prescribed velocity field will be considered here, so that: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 nullification, by two commands below, of high face porosities provides the independency of the slab sub-domains.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)

GROUP 11. Initialization of variable or porosity fields INIADD=F 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 U1AD=GRND {SCUF01} VELAD=-2.*U1[,,1] REGION(,NX-1) /IZ.EQ.2 V1AD=GRND {SCVF01} VELAD=-2.*V1[,,1] 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.

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 LFVELA =L0F(VELAD ) LFU1 =L0F(U1 ) DO 81001 IX=1 ,NX-1 IADD=NY*(IX-1) DO 81001 IY=1 ,NY I=IY+IADD 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 LFVELA =L0F(VELAD ) LFV1 =L0F(V1 ) DO 83001 IX=1 ,NX IADD=NY*(IX-1) DO 83001 IY=1 ,NY-1 I=IY+IADD L0VELA=LFVELA+I L0V1 =LFV1 +I 83001 F(L0VELA )=-2.*F(L0V1+NFM*(1-IZ)) ENDIF ENDIF ENDIF

Click here to return to "Contents"

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.

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.

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 IADD=NY*(IX-1) DO 96001 IY=IYF,IYL I=IY+IADD L0TEMP=LFTEMP+I L0VISL=LFVISL+I 96001 F(L0VISL )=1.0000E-03*EXP(-1.6000E+00*(F(L0TEMP)-0)) ENDIF ENDIF

Click here to return to "Contents"

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.

* 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)

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 IADD=NY*(IX-1) DO 98002 IY=IYF,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 99003 IY=IYF,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 10301 IY=IYF,IYL I=IY+IADD 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.

Click here to return to "Contents"

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

They cover:

- Flow field initialistions by parametric analytics
- Manipulating with initial fields
- MARKing sub-domains to create arbitrary initial fields
- Geometry initialisations

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.

Click here to return to "Contents"

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.

GROUP 13. Boundary conditions and special sources

REAL(RADCO,TMPX) RADCO=1.e-6; TMPX=100.;RG(1)=RADCO;RG(2)=TMPX PATCH(RAD,HIGH,1,NX,1,NY,NZ,NZ,1,1) {SORC01} CO=RG(1)*(RG(2)**2+H1**2)*(RG(2)+H1) {SORC01} VAL=RG(2) COVAL(RAD,h1,GRND,GRND)

c Source name: SORC01 IF(INDVAR.EQ.INAME('H1 ').AND.NPATCH.EQ.'RAD ') THEN LFCO =L0F(CO ) LFH1 =L0F(H1 ) DO 13701 IX=IXF,IXL IADD=NY*(IX-1) DO 13701 IY=IYF,IYL I=IY+IADD L0H1 =LFH1 +I 13701 F(LFCO +I)=RG(1)*(RG(2)**2+F(L0H1)**2)*(RG(2)+F(L0H1)) ENDIF IF(INDVAR.EQ.INAME('H1 ').AND.NPATCH.EQ.'RAD ') THEN LFVAL =L0F(VAL) DO 13801 IX=IXF,IXL IADD=NY*(IX-1) DO 13801 IY=IYF,IYL I=IY+IADD 13801 F(LFVAL+I)=RG(2) ENDIF

Click here to return to "Contents"

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.

{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.

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 IADD=NY*(IX-1) DO 19655 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 19656 IY=1 ,NY I=IY+IADD 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)

Click here to return to "Contents"

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:

- No relaxation for ISWEEP less and equal 100,
- Global relaxation for ISWEEP larger than 100 and less or equal 200 and
- 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.

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 IADD=NY*(IX-1) DO 13797 IY=IYF ,IYL I=IY+IADD L0TFAL=LFTFAL+I 13797 F(LFCO +I)=1./F(L0TFAL) ELSE LFCO =L0F(CO) DO IX=IXF,IXL IADD=NY*(IX-1) DO IY=IYF,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 13798 IY=IYF ,IYL I=IY+IADD L0TFAL=LFTFAL+I 13798 F(LFCO +I)=1./F(L0TFAL) ELSE LFCO =L0F(CO) DO IX=IXF,IXL IADD=NY*(IX-1) DO IY=IYF,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 13799 IY=IYF ,IYL I=IY+IADD L0TFAL=LFTFAL+I 13799 F(LFCO +I)=1./F(L0TFAL) ELSE LFCO =L0F(CO) DO IX=IXF,IXL IADD=NY*(IX-1) DO IY=IYF,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 19201 IY=2 ,3 I=IY+IADD 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 IADD=NY*(IX-1) DO 19202 IY=1 ,1 I=IY+IADD 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 IADD=NY*(IX-1) DO 19203 IY=1 ,1 I=IY+IADD 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 IADD=NY*(IX-1) DO 19204 IY=1 ,1 I=IY+IADD 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 IADD=NY*(IX-1) DO 19301 IY=1 ,NY I=IY+IADD 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

Click here to return to "Contents"

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.

{SC0301} RG(1)=AMAX1(RG(1),TAD) {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.

c Special calls name: SC0301 IF(ISTEP.GE.1 .AND.ISTEP.LE.LSTEP ) THEN IF(IZ.GE.1 .AND.IZ.LE.NZ ) THEN LFTAD =L0F(INAME('TAD ')) DO 19302 IX=1 ,NX IADD=NY*(IX-1) DO 19302 IY=1 ,NY I=IY+IADD L0TAD =LFTAD +I 19301 RG(1)=AMAX1(RG(1),F(L0TAD)) 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 IADD=NY*(IX-1) DO 19303 IY=1 ,NY I=IY+IADD 19302 VARMAX(INAME('TEMP'))=RG(1) ENDIF ENDIF

Click here to return to "Contents"

Click to return here after viewing pictures

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))

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 IADD=NY*(IX-1) DO 19601 IY=1 ,NY I=IY+IADD L0FS =LFFS +I L0TEMP =LFTEMP +I F(L0FS )=AMIN1(1.,AMAX1((1750-F(L0TEMP))/(1750-1310 1),0.0)) 19601 CONTINUE ENDIF ENDIF

Click here to return to "Contents"

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,

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.

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

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 IADD=NY*(IX-1) DO 90101 IY=IYF ,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 90102 IY=IYF ,IYL I=IY+IADD 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 IADD=NY*(IX-1) DO 19601 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 19602 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 19603 IY=1 ,NY I=IY+IADD 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 IADD=NY*(IX-1) DO 19604 IY=1 ,NY I=IY+IADD L0PROD=LFPROD+I L0OXID=LFOXID+I L0FUEL=LFFUEL+I 19604 F(L0PROD )=1.-F(L0OXID)-F(L0FUEL) ENDIF ENDIF