c
c filename gxpist.htm 070405
C**** SUBROUTINE GXPIST is called from GREX3, group 19, section 1,
C and is activated only when W1AD is set to ZMOVE and when AZW1
C is set to greater than zero in SATELLITE.
C Click here for the
C the PHOENICS Encyclopaedia description.
C
C.... The grid consists of two parts. The first part covers the zone
C from the top of the cylinder ( ie. z=0.0) to the top of the
C piston; this part contracts and expands as the piston moves
C up and down. The second part of the grid covers the bowl of
C the cylinder, & hence moves with the cylinder but does not
C expand or contract. The coding which follows assumes that the
C satellite-set grid distribution is the one appropriate for
C "bottom dead centre ( ie. TFIRST=0.0 ), even when a restart
C is in question, in which case TFIRST is not zero.
C
C.... The satellite-set parameters which enter here are:
C AZW1, the rotation rate in radians/sec;
C BZW1, the crank radius; and
C CZW1, the ratio of the connecting-rod length to the crank radius.
C IZW1, the iz location of the last slab in part one of an n-part
C moving grid.
C.... The library cases 327 and 330-331 make use of it.
C
SUBROUTINE GXPIST(ISTEP,NZ,TIM,DT,ISWEEP,LSWEEP,NPRINT,NTPRIN)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /LDATA/ LDAT(84)
LOGICAL LDAT,NULLPR
DIMENSION ZWO(2),ZWN(2),WAV(2),IZW(2)
SAVE IZW,DZBOWL,ZTDC,ZWO,ZWN,WAV
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
EQUIVALENCE (NULLPR,LDAT(32))
C
NAMSUB = 'GXPIST'
c.... precaution against incorrect nz setting
IF(NZ.LT.IZW1) RETURN
IF(GRNDNO(1,AZW1)) GO TO 100
IF(ISTEP.EQ.1) THEN
c.... initialization
IF(BZW1.LE.0.0.OR.CZW1.LE.0.0) THEN
CALL WRIT40('unsuitable settings for piston movement')
CALL SET_ERR(440,
* 'Error. Unsuitable settings for piston movement.',1)
C STOP
ENDIF
c.... the iz values of the ends of the first and second parts of the grid
IZW(1) = IZW1
IZW(2) = NZ
c
C.... ZWNZ is the EARTH index for distances of the HIGH faces from
C the z=0.0 plane.
c
C.... Function VARZ(INDVAR,IZZ) yields the F-array value, which related
C to z-wise variables, corresponding to INDVAR.
c.... the depth of the bowl
DZBOWL = VARZ(ZWNZ,NZ) - VARZ(ZWNZ,IZW1)
c
c.... the z value at "top-dead-centre", where part 1 of the grid has minimum
c depth
ZTDC = VARZ(ZWNZ,IZW1) - 2.0*BZW1
c
c.... the crank angle at the beginning of the first time step
ANGLE = AZW1 * (TIM-DT)
TERM = (SIN(ANGLE)/CZW1)**2
c
c.... the piston face position the
ZWO(1) = ZTDC + BZW1 * (1.+COS(ANGLE)+CZW1*(1.-SQRT(1.- TERM)))
c
c.... the bottom of bowl position then
ZWO(2) = ZWO(1) + DZBOWL
ENDIF
C.... conditions at the end of the current time step
ANGLE = AZW1*TIM
TERM = (SIN(ANGLE)/CZW1)**2
c
c.... new positions of piston face and bowl bottom
ZWN(1) = ZTDC + BZW1* (1.+COS(ANGLE)+CZW1*(1.-SQRT(1.- TERM)))
ZWN(2) = ZWN(1) + DZBOWL
c
c.... corresponding velocities
WAV(1) = (ZWN(1)-ZWO(1))/DT
WAV(2) = WAV(1)
GO TO 10
100 CONTINUE
c.................................. fixed piston velocity = bzw1
IF(ISTEP.EQ.1) THEN
IZW(1) = IZW1
IZW(2) = NZ
ZWN(2) = VARZ(ZWNZ,NZ)
ZWN(1) = VARZ(ZWNZ,IZW1)
WAV(1) = BZW1
WAV(2) = BZW1
ENDIF
ZWO(2) = ZWN(2)
ZWO(1) = ZWN(1)
ZWN(2) = ZWN(2) + BZW1*DT
ZWN(1) = ZWN(1) + BZW1*DT
10 CALL ZMOVE(IZW,WAV,ZWN,ZWO,2)
c.... print-out
IF(ISWEEP.EQ.LSWEEP .OR. MOD(ISWEEP,NPRINT).EQ.0) THEN
C.... ZWN is the grid location at the present time
C ZWO is the grid location at the previous time
C WAV is the velocity at which the grid moves
IF(.NOT.NULLPR) THEN
IF(MOD(ISTEP,NTPRIN).EQ.0) THEN
CALL WRIT40('two-part expanding-grid parameters at:')
CALL WRIT1R('time ',tim)
CALL WRIT2R('new zpis',zwn(1),'new zbwl',zwn(2))
CALL WRIT2R('old zbwl',zwo(1),'old zbwl',zwo(2))
CALL WRIT2R('vel pist',wav(1),'vel bowl',wav(2))
ENDIF
ENDIF
ENDIF
NAMSUB = 'gxpist'
END
c