Encyclopaedia Index

3. Examples of Use

  1. time-step specification
  2. expanding and contracting grid
  3. adaptive grid-size calculation
  4. analytical BFC grid generation
  5. convection-fluxes alteration
  6. a non-linear property correlation
  7. interphase transport
  8. devices for setting initial variables and porosity fields
  9. non-linear boundary conditions and source laws
  10. setting residual-reference values
  11. setting under-relaxation devices
  12. setting limits and increments
  13. output control
  14. Confined turbulent flame
  15. A shell-and-tube heat exchanger

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

 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

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


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.

This is what the generated grid will look like.

What the user puts into Q1

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.

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

        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"


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



   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.

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

        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"


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


Click here to return to "Contents"

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

    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.

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

       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"

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

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

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:



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

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

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)

What PLANT puts into GROUND



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"

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

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:

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

        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"


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



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

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

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

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

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

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:

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


Click here to return to "Contents"