SUBROUTINE INITAL C... C... M2 - UNSTEADY POISEUILLE FLOW IN A CIRCULAR TUBE C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NR=31) COMMON/T/ TAU, NSTOP, NORUN + /Y/ PHI(NR) + /F/ PHIT(NR) + /C/ DR, DRS, R(NR) + /I/ NG, IP C... C... SET NUMBER OF RADIAL GRID POINTS IF(NORUN.EQ.1)NG=11 IF(NORUN.EQ.2)NG=21 IF(NORUN.EQ.3)NG=31 C... C... SPATIAL GRID DR=1.0D0/DFLOAT(NG-1) DRS=DR**2 DO 1 I=1,NG R(I)=DFLOAT(I-1)*DR 1 CONTINUE C... C... INITIAL CONDITION DO 2 I=1,NG PHI(I)=0.0D0 2 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NR=31) COMMON/T/ TAU, NSTOP, NORUN + /Y/ PHI(NR) + /F/ PHIT(NR) + /C/ DR, DRS, R(NR) + /I/ NG, IP C... C... STEP THROUGH NG GRID POINTS DO 1 I=1,NG C... C... R = 0 IF(I.EQ.1)THEN PHIT(1)=4.0D0+4.0D0*(PHI(2)-PHI(1))/DRS C... C... R EQ 1 ELSE IF(I.EQ.NG)THEN PHI(NG)=0.0D0 PHIT(NG)=0.0D0 C... C... R NE 0 OR 1 ELSE PHIT(I)=4.0D0+(PHI(I+1)-2.0D0*PHI(I)+PHI(I-1))/DRS + +(1.0D0/R(I))*(PHI(I+1)-PHI(I-1))/(2.0D0*DR) END IF 1 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NR=31) COMMON/T/ TAU, NSTOP, NORUN + /Y/ PHI(NR) + /F/ PHIT(NR) + /C/ DR, DRS, R(NR) + /I/ NG, IP C... C... MONITOR SOLUTION WRITE(*,*)TAU,PHI(1) C... C... OPEN A FILE FOR MATLAB PLOTTING OF THE SOLUTION FOR THE 31 C... POINT GRID (NORUN = 31) IF((TAU.LT.0.001D0).AND.(NORUN.EQ.3))OPEN(1,FILE='m2.out') C... C... WRITE SOLUTION IF(NORUN.EQ.1)NGI=1 IF(NORUN.EQ.2)NGI=2 IF(NORUN.EQ.3)NGI=3 WRITE(NO,1)TAU,(PHI(I),I=1,NG,NGI) 1 FORMAT(/,' TAU = ',F6.2,/,' R', + 3X,' 0',3X,'0.1',3X,'0.2',3X,'0.3',3X,'0.4',3X,'0.5', + 3X,'0.6',3X,'0.7',3X,'0.8',3X,'0.9',3X,'1.0',/,' PHI', + 11F6.3,/) C... C... WRITE CENTERLINE VELOCITY TO MORE FIGURES FOR COMPARISON C... OF THE SOLUTIONS WRITE(NO,2)PHI(1) 2 FORMAT(' PHI(0) = ',F8.5,//) C... C... WRITE THE SOLUTION FOR MATLAB PLOTTING IP=IP+1 IF((NORUN.EQ.3).AND. + ((IP.EQ. 1).OR.(IP.EQ. 2).OR.(IP.EQ. 3).OR. + (IP.EQ. 4).OR.(IP.EQ. 5).OR.(IP.EQ. 7).OR. + (IP.EQ. 9).OR.(IP.EQ.11).OR.(IP.EQ.21)))THEN DO 3 I=1,NG WRITE(1,4)R(I),PHI(I) 4 FORMAT(F10.3,F10.5) 3 CONTINUE END IF RETURN END