driven by gravity.f90

伯肃叶流动用格子boltzmann方法去扩展新的研究方法

PROGRAM MAIN
PARAMETER(LY=11,LX=31)
INTEGER IS_SOLID_NODE(LY,LX)
REAL(8):: RHO(LY,LX),F(LY,LX,9),FTEMP(LY,LX,9),EX(9),EY(9),U_X(LY,LX),U_Y(LY,LX),X(LY,LX),Y(LY,LX),FEQ(LY,LX,9)
REAL(8):: TAU=1
REAL(8)::G=1.102E-3
REAL(8)::F1,F2,F3,RT0,RT1,RT2,UEQXIJ,UEQYIJ,UXSQ,UYSQ,UXUY5,UXUY6,UXUY7,UXUY8,USQ
!---SET SOLID NODES AT WALLS ON TOP AND BOTTOM---!
IS_SOLID_NODE=0
  DO I=1,LX
     IS_SOLID_NODE(1,I)=1
     IS_SOLID_NODE(LY,I)=1
  ENDDO
!---SET INITIAL DENSITY
 RHO=1.0
DO J=1,LY
   DO I=1,LX
      F(J,I,1)=(4./9.)*RHO(J,I)
      F(J,I,2)=(1./9.)*RHO(J,I)
      F(J,I,3)=(1./9.)*RHO(J,I)
      F(J,I,4)=(1./9.)*RHO(J,I)
      F(J,I,5)=(1./9.)*RHO(J,I)
      F(J,I,6)=(1./36.)*RHO(J,I)
      F(J,I,7)=(1./36.)*RHO(J,I)
      F(J,I,8)=(1./36.)*RHO(J,I)
      F(J,I,9)=(1./36.)*RHO(J,I)
   ENDDO
ENDDO
!---DEFINE LATTICE VELOCITY VECTORS---!
EX(1)=0
EY(1)=0
EX(2)=1
EY(2)=0
EX(3)=0
EY(3)=1
EX(4)=-1
EY(4)=0
EX(5)=0
EY(5)=-1
EX(6)=1
EY(6)=1
EX(7)=-1
EY(7)=1
EX(8)=-1
EY(8)=-1
EX(9)=1
EY(9)=-1

!--TIME LOOP--!
DO TS=1,3000
    WRITE(*,*) TS
!---COMPUTATING MACROSCOPIC DENSITY,RHO,AND VELOCITY,U=U(UX,UY)---!
DO J=1,LY
   DO I=1,LX
      U_X(J,I)=0.0
	  U_Y(J,I)=0.0
	  RHO(J,I)=0.0
        IF(IS_SOLID_NODE(J,I).EQ.0)THEN
           DO K=1,9
	          RHO(J,I)=RHO(J,I)+F(J,I,K)
		      U_X(J,I)=U_X(J,I)+EX(K)*F(J,I,K)
		      U_Y(J,I)=U_Y(J,I)+EY(K)*F(J,I,K)
           ENDDO
	           U_X(J,I)=U_X(J,I)/RHO(J,I)
	           U_Y(J,I)=U_Y(J,I)/RHO(J,I)
	    ENDIF
	!---ADD SPACE MATRICIES FOR PLOTTING---!
	  X(J,I)=I
	  Y(J,I)=J
    ENDDO
ENDDO

!---COMPUTE THE EQUILIBRIUM DISTRIBUTION FUNCTION,FEQ---!
F1=3.
F2=9./2.
F3=3./2.
DO J=1,LY
   DO I=1,LX
        IF(IS_SOLID_NODE(J,I).EQ.0)THEN

		   RT0=(4./9.)*RHO(J,I)
		   RT1=(1./9.)*RHO(J,I)
		   RT2=(1./36.)*RHO(J,I)

		   UEQXIJ=U_X(J,I)+TAU*G
		   UEQYIJ=U_Y(J,I)

		   UXSQ=UEQXIJ*UEQXIJ
		   UYSQ=UEQYIJ*UEQYIJ

		   UXUY5= UEQXIJ+UEQYIJ
		   UXUY6=-UEQXIJ+UEQYIJ
		   UXUY7=-UEQXIJ-UEQYIJ
		   UXUY8= UEQXIJ-UEQYIJ
          
		   USQ=UXSQ+UYSQ

		   FEQ(J,I,1)=RT0*(1-F3*USQ)

		   FEQ(J,I,2)=RT1*(1+F1*UEQXIJ+F2*UXSQ-F3*USQ)
		   FEQ(J,I,3)=RT1*(1+F1*UEQYIJ+F2*UYSQ-F3*USQ)
		   FEQ(J,I,4)=RT1*(1-F1*UEQXIJ+F2*UXSQ-F3*USQ)
		   FEQ(J,I,5)=RT1*(1-F1*UEQYIJ+F2*UYSQ-F3*USQ)

		   FEQ(J,I,6)=RT2*(1+F1*UXUY5+F2*UXUY5*UXUY5-F3*USQ)
		   FEQ(J,I,7)=RT2*(1+F1*UXUY6+F2*UXUY6*UXUY6-F3*USQ)
		   FEQ(J,I,8)=RT2*(1+F1*UXUY7+F2*UXUY7*UXUY7-F3*USQ)
		   FEQ(J,I,9)=RT2*(1+F1*UXUY8+F2*UXUY8*UXUY8-F3*USQ)
		ENDIF
	ENDDO
ENDDO


!---COLLISION STEP---!
DO J=1,LY
   DO I=1,LX
      
	  IF(IS_SOLID_NODE(J,I).EQ.1)THEN
	    !---STANDARD BOUNCEBACK---!
		 TEMP=F(J,I,2)
		 F(J,I,2)=F(J,I,4)
         F(J,I,4)=TEMP

		 TEMP=F(J,I,3)
		 F(J,I,3)=F(J,I,5)
         F(J,I,5)=TEMP

		 TEMP=F(J,I,6)
		 F(J,I,6)=F(J,I,8)
         F(J,I,8)=TEMP

		 TEMP=F(J,I,7)
		 F(J,I,7)=F(J,I,9)
         F(J,I,7)=TEMP
	    ELSE
	      !---REGULAR COLLISION---!
		    DO K=1,9
		       F(J,I,K)=F(J,I,K)-(F(J,I,K)-FEQ(J,I,K))/TAU
			ENDDO
	    ENDIF
	ENDDO
ENDDO

!---STREAMING STEP:SUBTLE CHANGES TO PERIODICITY HERE DUE TO INDEXING---!
   DO J=1,LY
      IF(J.GT.1)THEN
	    JN=J-1
	  ELSE
	    JN=LY
	  ENDIF
    
	 IF(J.LT.LY)THEN
	    JP=J+1
	  ELSE
	    JP=1
	  ENDIF

	DO I=1,LX
	   IF(I.GT.1)THEN
	    IN=I-1
	  ELSE
	    IN=LX
	  ENDIF

	  IF(I.LT.LX)THEN
	    IP=I+1
	  ELSE
	    IP=1
	  ENDIF

	  FTEMP(J,I,1)=F(J,I,1)
	  FTEMP(J,IP,2)=F(J,I,2)
	  FTEMP(JP,I,3)=F(J,I,3)
	  FTEMP(J,IN,4)=F(J,I,4)
	  FTEMP(JN,I,5)=F(J,I,5)
	  FTEMP(JP,IP,6)=F(J,I,6)
	  FTEMP(JP,IN,7)=F(J,I,7)
	  FTEMP(JN,IN,8)=F(J,I,8)
	  FTEMP(JN,IP,9)=F(J,I,9)
	ENDDO
ENDDO

   F=FTEMP
!---END TIME LOOP---!
    ENDDO   !---与 TIME LOOP 处的DO相呼应---!

	OPEN(UNIT=20,FILE='VELOCITY.DAT')
	 WRITE(20,*) 'VARIABLES = X, Y, U, V'                                                   
    WRITE(20,*) 'ZONE I=',LX,',J=',LY,' F=point'
     
	DO J=1,LY
	DO I=1,LX

	   WRITE(20,"(4F12.6)"),X(J,I),Y(J,I),U_X(J,I),U_Y(J,I)
    ENDDO
	ENDDO


    


	ENDPROGRAM