PARAMETER (MDX=201) COMMON /SETUPS/ G,R1,P1,E1,U1,R4,P4,E4,U4,XD0 COMMON /CMPCND/ MX,DX,X(MDX),CFL,DT,NLAST,DELTA COMMON /FLWDAT/ R(MDX),P(MDX),E(MDX),U(MDX),Q(3,MDX),DQ(3,MDX) DIMENSION UM(MDX),CM(MDX),HM(MDX) DIMENSION AL1(MDX),AL2(MDX),AL3(MDX) DIMENSION PH1(MDX),PH2(MDX),PH3(MDX) DIMENSION TV1(MDX),TV2(MDX),TV3(MDX) WRITE(6,*) '>> SOLVE SHOCK-TUBE PROBLEM USING TVD METHOD' C......(1) SHOCK-TUBE SETUP G = 1.4 R1 = 0.25 P1 = 0.2154 E1 = P1/((G-1.)*R1) U1 = 0.0 R4 = 1.0 P4 = 1.0 E4 = P4/((G-1.)*R4) U4 = 0.0 XD0 = 0.0 X40 = -0.6 X10 = 0.6 WRITE(6,*) ' ' WRITE(6,*) '>> SHOCK-TUBE SETUP' WRITE(6,*) ' HIGH PRESSURE SIDE : PRESSURE = ', P4 WRITE(6,*) ' DENSITY = ', R4 WRITE(6,*) ' LOW PRESSURE SIDE : PRESSURE = ', P1 WRITE(6,*) ' DENSITY = ', P4 C......(2) GRID GENERATION MX = 121 DX = (X10-X40)/FLOAT(MX-1) DX2 = DX * 2. DXH = DX / 2. DO 200 I=1,MX+1 X(I)=X40+DX*FLOAT(I-1) 200 CONTINUE C......(3) TIME SETUP SIZE WRITE(6,*) '>> INPUT CFL NUMBER ===>' READ(5,*) CFL VREF=1.0 DT = CFL*DX/VREF WRITE(6,*) 'TIME STEP SIZE =', DT WRITE(6,*) '>> INPUT NO. OF TIME STEPS ===>' READ(5,*) NLAST C......(4) TVD SOLUTION C......4-1 INITIAL CONDITION WRITE(6,*) ' ' WRITE(6,*) '>> NUMERICAL SIMULATION STARTS !!' DO 410 I=1,MX IF(X(I).LT.XDO) THEN R(I) = R4 P(I) = P4 E(I) = E4 U(I) = U4 ELSE R(I) = R1 P(I) = P1 E(I) = E1 U(I) = U1 END IF Q(1,I)=R(I) Q(2,I)=R(I)*U(I) Q(3,I)=R(I)*E(I)+0.5*R(I)*U(I)**2 410 CONTINUE C......4-2 TIME MARCHING DO 4200 N=1,NLAST TIME=DT*FLOAT(N) IF((N/10)*10.EQ.N) WRITE(6,*) ' ... STEP = ',N C <1> I+1/2 VALUES DO 4210 I=1,MX-1 C <== ROE'S AVERAGING DASH=SQRT(R(I+1)/R(I)) UM(I)=(DASH*U(I+1)+U(I))/(DASH+1.) HM(I)=(DASH*(Q(3,I+1)+P(I+1))/R(I+1) & + (Q(3,I)+P(I))/R(I))/(DASH+1.) CM2=(G-1.)*(HM(I)-0.5*UM(I)**2) CM(I)=SQRT( & AMAX1(CM2, AMIN1(G*(G-1.)*E(I), G*(G-1.)*E(I+1)))) DQ1=Q(1,I+1)-Q(1,I) DQ2=Q(2,I+1)-Q(2,I) DQ3=Q(3,I+1)-Q(3,I) C <== COMPUTR R^-1 B2=(G-1.)/CM(I)**2 B1=B2*UM(I)**2/2. A11=0.5*(B1+UM(I)/CM(I)) A12=0.5*(-B2*UM(I)-1./CM(I)) A13=0.5*B2 A21=1.-B1 A22=B2*UM(I) A23=-B2 A31=0.5*(B1-UM(I)/CM(I)) A32=0.5*(-B2*UM(I)+1./CM(I)) A33=0.5*B2 C <== COMPUTE A=R^-1*DQ AL1(I)=A11*DQ1+A12*DQ2+A13*DQ3 AL2(I)=A21*DQ1+A22*DQ2+A23*DQ3 AL3(I)=A31*DQ1+A32*DQ2+A33*DQ3 4210 CONTINUE C <2> LIMITER FUNCTION C MINMOD(X,Y,Z)=S*MAX[0,MIN(S*X,S*Y,S*Z)], S=SGN(X) C <== ENTROPY CORRECTION PARAMETER DELTA=0.15 DO 4220 I=2,MX-2 RAMDA=UM(I)-CM(I) FF1=(DT/DX)*RAMDA**2 FF2=ABS(RAMDA) IF(FF2.LT.DELTA) FF2=0.5*(RAMDA**2+DELTA**2)/DELTA S=SIGN(1. , AL1(1-I)) QQ=S*AMAX1(0., AMIN1(S*2.*AL1(I-1), S*2.*AL1(I), & S*2.*AL1(I+1), S*0.5*(AL1(I-1)+AL1(I+1)))) PH1(I)=-FF1*QQ-FF2*(AL1(I)-QQ) RAMDA=UM(I) FF1=(DT/DX)*RAMDA**2 FF2=ABS(RAMDA) IF(FF2.LT.DELTA) FF2=0.5*(RAMDA**2+DELTA**2)/DELTA S=SIGN(1. , AL2(I-1)) QQ=S*AMAX1(0., AMIN1(S*2.*AL2(I-1), S*2.*AL2(I), & S*2.*AL2(I+1), S*0.5*(AL2(I-1)+AL2(I+1)))) PH2(I)=-FF1*QQ-FF2*(AL2(I)-QQ) RAMDA=UM(I)+CM(I) FF1=(DT/DX)*RAMDA**2 FF2=ABS(RAMDA) IF(FF2.LT.DELTA) FF2=0.5*(RAMDA**2+DELTA**2)/DELTA S=SIGN(1. , AL3(I-1)) QQ=S*AMAX1(0., AMIN1(S*2.*AL3(I-1), S*2.*AL3(I), & S*2.*AL3(I+1), S*0.5*(AL3(I-1)+AL2(I+1)))) PH3(I)=-FF1*QQ-FF2*(AL3(I)-QQ) 4220 CONTINUE C <3> NUMERICAL FLUX DO 4230 I=2,MX-2 TVD1=PH1(I)+PH2(I)+PH3(I) TVD2=UM(I)*(PH1(I)+PH2(I)+PH3(I))-CM(I)*(PH1(I)-PH3(I)) TVD3=(HM(I)-UM(I)*CM(I))*PH1(I)+0.5*UM(I)**2*PH2(I) & +(HM(I)+UM(I)*CM(I))*PH3(I) TV1(I)=(Q(2,I)+Q(2,I+1)+TVD1)/2. TV2(I)=(Q(2,I)*U(I)+P(I)+Q(2,i+1)*U(I+1)+P(I+1)+TVD2)/2. TV3(I)=((Q(3,I)+P(I))*U(I)+(Q(3,I+1)+P(I+1))*U(I+1) & +TVD3)/2. 4230 CONTINUE C <4> DIFFERENCING DO 4240 I=3,MX-2 DQ(1,I)=-(DT/DX)*(TV1(I)-TV1(I-1)) DQ(2,I)=-(DT/DX)*(TV2(I)-TV2(I-1)) DQ(3,I)=-(DT/DX)*(TV3(I)-TV3(I-1)) 4240 CONTINUE C <5> UP-DATE DO 4250 I=3, MX-2 Q(1,I)=Q(1,I)+DQ(1,I) Q(2,I)=Q(2,I)+DQ(2,I) Q(3,I)=Q(3,I)+DQ(3,I) 4250 CONTINUE C <6> SET BOUNDARY CONDITIONS Q(1,1)=R4 Q(2,1)=R4*U4 Q(3,1)=R4*E4+0.5*R4*U4**2 Q(1,2)=R4 Q(2,2)=R4*U4 Q(3,2)=R4*E4+0.5*R4*U4**2 Q(1,MX-1)=R1 Q(2,MX-1)=R1*U1 Q(3,MX-1)=R1*E1+0.5*R1*U1**2 Q(1,MX)=R1 Q(2,MX)=R1*U1 Q(3,MX)=R1*E1+0.5*R1*U1**2 C <7>PRIMITIVE VARIABLES DO 4270 I=1,MX R(I)=Q(1,I) U(I)=Q(2,I)/R(I) E(I)=Q(3,I)/R(I)-0.5*U(I)**2 P(I)=(G-1.)*R(I)*E(I) 4270 CONTINUE C 4200 CONTINUE C......(5) WRITE TVD SOLUTION IN FILE WRITE(6,*) ' ' WRITE(6,*) '>> WRITE CFD SOLUTION ON FILE.' OPEN(UNIT=60,FILE='SOCK.DATA',FORM='FORMATTED') WRITE(60,*) MX,CFL,TIME,NLAST WRITE(60,*) ' x Density Pressure e Velocity' DO 500 I=1,MX WRITE(60,501) X(I),R(I),P(I),E(I),U(I) 500 CONTINUE 501 FORMAT(5F15.4) CLOSE(UNIT=60) C WRITE(6,*) ' ' WRITE(6,*) '>> COMPLETED.' C C STOP END