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