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 '/phoenics/d_includ/farray' INCLUDE '/phoenics/d_includ/grdloc' INCLUDE '/phoenics/d_includ/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