2 С SHOCK COMPUTES THE PROPAGATION OF A SHOCK USING

3 С HACCORMACK OR LAX-WENDROFF SCHEMES

4 С WITH ARTIFICIAL VISCOSITY OR FCT SMOOTHING

6 DIMENSION X(lOl),P(101>.RH(lOl),U(101).T(lOl),UEX(101)

7 DIMENSION Q(101,3).QD{101.3).r(10i,3),DF(101,3),DL(3),DbM(3)

8 OPEN (1,FILE= SHOCK. DAT*)

9 OPEN(б,FILE= SHOCK.OUT)

10 READ(1,1)IME,IFCT,I??.,NX,NT,DT,GAM,PRAT

11 READ(1,2)ETA,ET1,E72,ENL,SHST

12 1 FO?HAT(5I5,3F5.3)

13 2 FORMAT{3F8.5,2F5.3)

14 С

15 ?I = 3.1415927

16 NXM = N7. - 1

17 ANX = NXM

18 DX = l./AMX

19 JSST = SHST/DX .20 DTR = DT/DX

21 CMM = GAM - 1.

22 GMP s GAM + 1.

23 GMR = CMM/GMP

24 SHS = 0.5MGMM + GMP*PRAT)/GAM

25 SHS = SQRT(SHS)

26 С

27 IFdME .EQ. 2)WRITE{6,4)

28 IFdME .EQ. 1)WRITE{6,3)

29 IFCIFCT.EQ. 1)WRITE(6,5)ETA.ET1,ET2

30 WRITE{6,6)NX,NT,DX,DT,ENL

31 WRITE(6.7)CAM,PRAT,SHS,SHST

32 3 FORMATC ONE-DIMENSIONAL SHOCK PROPAGATION, MACCORMACK SCHEME)

33 4 FORMATC ONE-DIMENSIONAL SHOCK PROPAGATION, LAX-WENDROFF*

34 1, SCHEME)

35 5 FORMATC FLUX CORRECTED TRANSPORT, ETA,ET1,ET2=,3F8.5)

36 6 FORMATC NX=M3. NT=M3, DX=,F5.3, DT=,F5.3, ENL=,P5.3)

37 7 FORMATC GAM=,F5.3, Pl/P2=,F5.2, SH/SP=,F5.3, SHST=,

38 lr5.3,/}

39 С

40 С SET i:UTIAL CONDITIONS

41 С

42 DUM = GMR + PRAT

43 RHD = DUM/(1. + GMR*PRAT)

44 DIM = SQRT{2.*GAM/GMP/DUM)

45 UD = (PRAT - l.)*DIH/GAM

46 DO 8 J = l.JSST

47 U(J) = UD

48 RH(J) = RHD

49 P(J) = PRAT

50 T(J) P(J)/RH(J)

51 8 CONTINUE

52 JHR = JSST + 1

53 DO 9 J = JHR,NX

54 U(J) = 0.

55 RH(J) = 1.

56 P(J) = 1.

57 T(J) = 1.

58 9 CONTINUE

59 TIM = 0.

60 N 0

61 С

62 С SET INITIAL Q AND F

63 С

€4 БО 10 J = 1,NX

65 AJ J - 1

66 X(J) = AJ*DX

67 Q(J,1) RH(J)

68 Q(J,2) RH(J)*U(J)

69 Q(J,3) = P(J)/GAM/GMH + 0.5*RH(J)*U(J)*U(J)

70 F(J,1) =.Q(J,2)

71 r(J,2) = P(J)/GAM + U(J)*Q(J,2)

72 F(J,3) = (P(J)/GHH + 0.5*U(J)*Q(J,2))*U(J)

73 10 CONTINUE

74 WRITE(6,21)N,TIM

75 VRITE(6,22)(X(J),a=l,NX)

76 WRITE(6,23)(U(J),J=1.NX)

77 WRITE(6,24)(RH(J),J 1,NX)

78 ¥RITE(6,25)(P(J),J=1,NX)

79 WRITE(6,26)(T(J),a=l,NX)

80 С

81 С ADVANCE SOLUTIOM IH TIME

83 DO 30 N = 1,NT

84 AN N

85 IF (IFCT .NE. DGOTO 12

86 ENU = ETA + 0.25*ET1*((U(1)+U(2))*DTR)**2

87 DO 11 К = 1,3

88 11 DF(bK) = ENU*(Q(2,K)-Q(1,K))

89 12 CONTINUE

90 С

91 С OBTAIN HALF-STEP SOLUTION

92 С -

93 DO 18 J = 2,NXM

94 p. DO 13 К 1,3

95. IFdME .EQ. 1)QD(J,K) = Q(J,K) - ©TR* (Г (J+1,K) - F(J,K))

96 IFdHE .EQ. 2)QD(J,K) 0.5*(Q(J,K)+Q(J+1,K)) - 0.50*DTR*

97 1(F(J+1,K) - F(J,K))

98 IFdFCT .EQ. DENU = ETA+0.25*ETl* ( (UTJ)+U(J+1))*DTR)**2

99 1 IFdFCT .EQ, 1)DF(J,K) = ENU*(Q(J+1,K) - Q(J,K))

100 iLl3 CONTINUE

101 F(J,1) = QD(J,2)

102 UD = QD(J,2)/QD(J,1)

103 PD (QD(J,3) - 0.5*UD*QD(J,2))*GAM*GHM

104 F(J,2) = PD/GAM + UD*QD(J,2)

105 F(J,3) (PD/GMM + 0.5*UD*QD(J,2))*UD

106 С

107 С OBTAIN FULL-STEP SOLUTION

108 С

109 DO 14 К = 1,3

110 IFdME .EQ. 1)Q(J,K) = 0.5*(Q(J,K)+QD(J,K)) - 0.5*DTR*(F(JД)

111 1 - F(J-1,K))

112 IFdME .EQ. 2)Q(J,K) = Q(J,K) - DTR*(F(J,K) - F(J-1,K))

113 14 CONTINUE

114 С

115 С ARTIFICAL VISCOSITY SMOOTHING

116 С

117 IFdFCT .EQ. DGOTO 18

118 IF(J .NE. 2}G0T0 16

119 DUC ABS(U{2) - U(l))

120 IF(DUC .LT. 1.0E-06)DUC = l.OE-06

121 DO 15 К = 1,3

122 DLM(K) = Q(2,K) - Q(1,K)

123 IF(ABS(DLM(K)) .LT. 1.OE-06)DLM(K) - 1.0E-06*SIGN(1.0 DLM(K))

124 DUC = ABS(DLM(K))

125 15 DLM(K) = DUC*DLM(K)

126 16 DUC = ABS(U(J+1)-U(J))

127 IF(DUC .LT. l.OE-06)DUC = l.OE-06

Рис. 14.17 (продолжение).

DO 17 К = 1,3

DLCK) = Q(J+1,K) - Q(J,K)

IF(ABS(DL(K)) .LT.1.0E-06)DL(K)=1.0E-06*SIGN(1.0 DL(K)J DUC ABS(DL(K}) DLCK) = DUC*DL{K)

Q(J,K) = Q{J.K) + ENL*DTR*(DL(K) - DLH(K)) DLHCK) = DL(K)

FCT SMOOTHING OBTAIN RH,U,P AND F

IFdFCT .EQ. DCALL FCT(NXM,ETA,ET2,DTR,Q,DF,U)

DO 19 J = 2,NXM

RH(J) = Q(J,1)

U{J) = Q{J.2)/Q(J,1)

P(J) = (Q(J,3) - 0.5*U(J)*Q(J,2)) GAM*GHH

T(J) = P(J)/RH(J)

F(J,1) = RH(J)*U(J)

F{J,2) = P(J)/GAM + RH(J)*U(J)*U(J)

F(J,3) = (P(J)/GMM + 0.5*RH(J)*U(J)*U(J))*U(a)

TIM = AN*DT

IF(N .EQ. NT)GOTO 20

IF(IPR .EQ. 0)GOTO 30

WRITE(6,21)N,TIM

FORMAT(/, N=M3, TIM ,E10,3) WRITE(6,22)(X(J),J 1,NX) WRITE(6,23)(U(J).J=1,NX) VRITE(6,24) (RH(J), J=1,NX) VRITE(6,25)(P(J),J 1,NX) WRITE(6.26)(T(J),J=1,NX) FORMATC X M2F6.3) FORMATC U=M2F6.3) FORMATC RH=M2F6.3) FORMATC P=M2F6.3) FORMATC T=M2F6.3) IFdPR .LE. DGOTO 30 DO 27 К = 1,3

WRITE(6,28)(Q(J,K),J=1,NX) WRITE(6,29)(F(J,K),J=1,NX) CONTINUE

FORMATC Q=M2F6.3) FORMATC F=M2F6.3)

EXACT SOLUTION, UEX

SHFN = SHST + SHS*TIM

JSST = SHFN/DX +1.0

DO 31 J = 1,JSST

IF(J .LE. JSST)UEX(J) = U(l)

IF(J .GT. JSST)UEX(J) = U(NX)

CONTINUE

VRITE(6,32)(UEX(J),J=1,NX) FORMATC UEX=M2F6.3) STOP END