shock_tube.f90

一维激波管程序

      PROGRAM RIEMANN
	  IMPLICIT NONE
	  REAL*8,PARAMETER :: T=0.3,DX=0.002,GAMA=1.4,CFL=0.5
	  INTEGER,PARAMETER :: SN=500,TMAX=900
	  REAL*8  U(-SN:SN,0:TMAX,3),U1(-SN:SN,0:TMAX,3),&
	          V(-SN:SN,0:TMAX),P(-SN:SN,0:TMAX),F(3),F1(3),&
			  AV(3),A(-SN:SN),LAMDA,Q,DT,ST,VMAX,SS
	  INTEGER I,J,TN
	  !u(i,j,1)--roul(i,j),u(i,j,2)--(roul*v)(i,j)
	  !u(i,j,3)--e=p/(gama-1)+1./2*roul*u**2
	  !AV--artificial viscidity, a--sound speed
	  !VMAX--the max speed of the shock,F--flux,LAMDA--dt/dx

      !INITIAL VALUE
	  VMAX=0
	  DO I=-SN,SN
	     IF(I<=0) THEN
		    V(I,0)=0.
			P(I,0)=1.
			U(I,0,:)=(/1.,0.,2.5/)
		  ELSE
		    V(I,0)=0.
			P(I,0)=0.1
			U(I,0,:)=(/0.125,0.,0.25/)
		  ENDIF
		  U1(I,0,:)=U(I,0,:)
		  A(I)=SQRT(GAMA*P(I,0)/U(I,0,1))
		  VMAX=MAX(VMAX,ABS(V(I,0)+A(I)),ABS(V(I,0)-A(I)))
	  END DO

	!CACULATE
	  DT=CFL*DX/VMAX
	  TN=NINT(T/DT)
	  LAMDA=DT/DX
	  !WRITE THE NUMERICAL INITIAL VALUE "dt,tn,lamda" into file
	  OPEN(UNIT=15,FILE='INI_VALUE.DAT')
	      WRITE(15,*) 'FIRST VALUE:'
	      WRITE(15,*)
	      WRITE(15,*) 'DT=',DT
	      WRITE(15,*) 'TN=',TN
  	      WRITE(15,*) 'LAMDA=',LAMDA
	  CLOSE(UNIT=15)
	  OPEN(UNIT=16,FILE='FINAL_VALUE.DAT')
	  !main process
	  ST=0.           !the sum time
      ss=vmax*dt      !the sum space of the shock propagating
	  DO I=0,TMAX-1
	     VMAX=0.
         ST=ST+DT
         CALL FLUX(U(-SN,I,:),U(-SN+1,I,:),F,F1)
		 U1(-SN,I+1,:)=U(-SN,I,:)-LAMDA*(F1-F)
		 DO J=-SN+1,SN-1
		    CALL FLUX(U(J,I,:),U(J+1,I,:),F,F1)
			Q=0.5*ABS((P(J+1,I)-2*P(J,I)+P(J-1,I))/(P(J+1,I)+2*P(J,I)+P(J-1,I)))
		!between the space of the shock propagating add the artificial viscidity 
			    AV(:)=Q*(U(J+1,I,:)-2*U(J,I,:)+U(J-1,I,:))
			U1(J,I+1,:)=U(J,I,:)-LAMDA*(F1-F)+AV
			CALL FLUX(U1(J-1,I+1,:),U1(J,I+1,:),F,F1)
		    	AV=Q*(U1(J+1,I,:)-2*U1(J,I,:)+U1(J-1,I,:))
			U(J,I+1,:)=1./2*(U(J,I,:)+U1(J,I+1,:))-1./2*LAMDA*(F1-F)+AV
			V(J,I+1)=U(J,I+1,2)/U(J,I+1,1)
			P(J,I+1)=(GAMA-1)*(U(J,I+1,3)-1./2*U(J,I+1,2)*V(J,I+1))
			A(J)=SQRT(GAMA*P(J,I+1)/U(J,I+1,1))
			VMAX=MAX(VMAX,ABS(V(J,I+1)+A(J)),ABS(V(J,I+1)-A(J)))
		 END DO
		 U(-SN,I+1,:)=U(-SN+1,I+1,:)
		 V(-sn,I+1)=U(-sn,I+1,2)/U(-sn,I+1,1)
		 P(-sn,I+1)=(GAMA-1)*(U(-sn,I+1,3)-1./2*U(-sn,I+1,2)*V(-sn,I+1))
		 U(SN,I+1,:)=U(SN-1,I+1,:)
		 V(sn,I+1)=U(sn,I+1,2)/U(sn,I+1,1)
		 P(sn,I+1)=(GAMA-1)*(U(sn,I+1,3)-1./2*U(sn,I+1,2)*V(sn,I+1))
		 DT=CFL*DX/VMAX
		 LAMDA=DT/DX
		 ss=ss+vmax*dt
		 IF(ST>=T) THEN
	        WRITE(16,*) 'FINAL VALUE:'
		    WRITE(16,*) 'FINAL TIME:', ST
			EXIT
		 ENDIF
      END DO
	  tn=i
	  WRITE(16,*) 'ST=',ST
	  WRITE(16,*) 'TN=',TN
	  WRITE(16,*) 'DT=',DT
	  WRITE(16,*) 'LAMDA=',LAMDA
	 !WRITE DATA
	    OPEN(UNIT=10,FILE='data.dat')
		   write(10,*)"variables=x,p,roul,v"
		   write(10,*)"zone i=",2*sn
		   DO J=-SN,SN
	          write(10,*) j*dx,p(J,TN),U(J,TN,1),v(J,TN)
    	   END DO
		CLOSE(UNIT=10)
        !OPEN(UNIT=11,FILE='VELOCITY.dat')
		 !  DO J=-SN,SN
		     ! write(11,*) j*dx,v(j,tn)
		  ! END DO
	!	CLOSE(UNIT=11)
	!	OPEN(UNIT=12,FILE='PRESSURE.DAT')
		!   WRITE(12,*) "PRESSURE:"
		!   DO J=-SN,SN
		 !     WRITE(12,*) J*DX, P(J,TN)
		 !  END DO
        !CLOSE(UNIT=12)
	  END
    
	!SUBROUTINE
	   !caculate the flux  
       subroutine FLUX(u0,u1,f,f1)
	   IMPLICIT NONE	
	   real*8 u0(3),u1(3),f(3),f1(3)
	   real*8 vm,b(3,3)
	   
	   vm=u0(2)/u0(1)
	   call trans_matrix(vm,b)
	   call multi_matrix(u0,b,f)
	   vm=u1(2)/u1(1)
	   call trans_matrix(vm,b)
	   call multi_matrix(u1,b,f1)
	   return
	   end

	!subroutine 2
	   subroutine trans_matrix(vm,b)
	   implicit none
	   real*8 vm,b(3,3)
	   real,parameter :: gam=1.4
	   b=0
	   b(1,1)=vm
	   b(2,2)=(3-gam)/2*vm
	   b(2,3)=gam-1
	   b(3,2)=(1-gam)/2*vm**2
	   b(3,3)=gam*vm
	   return
	   end

	!subroutine 3
	   subroutine multi_matrix(u,b,f)
	   implicit none
	   real*8 u(3),b(3,3),f(3)
	   integer i,j
	   f=0
	   do j=1,3
	      do i=1,3
		     f(j)=f(j)+b(j,i)*u(i)
		  end do
		end do
		return
		end