qbalance.for

关于二维射流计算源程序

	SUBROUTINE      QBALANCE
	 USE VOF2D
      IMPLICIT        REAL*8(A-H,O-Z)
      CHARACTER*15     FILENAME3,FILENAME4
!	INCLUDE         'vof.INC'
      real*8     sumu(6,mj), sumvis(6,mj), sumTK(6,mj), sumTC(6,mj)
      real*8     FAI(MI,MJ)
	real*8     EVEETA(MI),evebed(mi)
c----------------对流入和流出量进行累加,并且求差值	  
          
c----------在0时刻给变量赋值	
      IF(ISTYLE.ge.2)THEN
        IF(T.EQ.0)THEN  
		  DO I=1,IMAX
	       DO J=1,JMAX
	         EVEU(I,J)=0
	         EVEV(I,J)=0
	         EVETC(I,J)=0.0
	         EVETC1(I,J)=0.0
	         EVETC2(I,J)=0.0
	         EVEUN(I,J)=U(I,J)
	         EVEVN(I,J)=V(I,J)
	         EVETCN(I,J)=TC(I,J)
	         EVETCN1(I,J)=TC1(I,J)
	         EVETCN2(I,J)=TC2(I,J)
	       ENDDO
	       EVEETA(I)=0.
             EVEBED(I)=0.
	      ENDDO
	       

	    DO I=1,NNETA
		 DO J=2,JMAX
		   SUMU(I,J)=0.
		   SUMVIS(I,J)=0.
	       SUMTK(I,J)=0.
	       SUMTC(I,J)=0.
		ENDDO

		ENDDO
	  ELSE

	DO I=1,IM1
	DO J=1,JM1
	EVEU(I,J)=EVEU(I,J)+U(I,J)*F(I,J)
	EVEV(I,J)=EVEV(I,J)+V(I,J)*F(I,J)
	EVETC(I,J)=EVETC(I,J)+TC(I,J)*F(I,J)
	EVETC1(I,J)=EVETC1(I,J)+TC1(I,J)*F(I,J)
	EVETC2(I,J)=EVETC2(I,J)+TC2(I,J)*F(I,J)
	ENDDO
	EVEETA(I)=EVEETA(I)+ETAA(I)
	EVEBED(I)=EVEBED(I)+USBED(I)
	ENDDO
		  
		DO I=1,NNETA
		  II=IETA(I)
	    
		 DO J=2,JMAX
		  IF(NF(II,J).NE.6)THEN
		   SUMU(I,J)=SUMU(I,J)+U(II,J)*F(II,J)
		   SUMVIS(I,J)=SUMVIS(I,J)+VIS(II,J)*F(II,J)
	       SUMTK(I,J)=SUMTK(I,J)+TK(II,J)*F(II,J)
	       SUMTC(I,J)=SUMTC(I,J)+TC(II,J)*F(II,J)
	      ENDIF
		 ENDDO
		ENDDO
  
		DO N1=1,NSUM
          IF(T.GT.(TWAVE*(N1+0.5)-0.5*DELT).
     *		 AND.T.LT.(TWAVE*(N1+0.5)+0.5*DELT))THEN	
     	            
		  CALL CALTRAGICY(FAI,EVEU,EVEV)      
 

		 
		   DO I=2,IM1
	        DO J=2,JM1
		      EVEU(I,J)=EVEU(I,J)/NTEVE
		      EVEV(I,J)=EVEV(I,J)/NTEVE
	          EVETC(I,J)=EVETC(I,J)/NTEVE
	          EVETC1(I,J)=EVETC1(I,J)/NTEVE
	          EVETC2(I,J)=EVETC2(I,J)/NTEVE
!	          WRITE(61,121)(X(I)-X(IETA(2)))/FLHT0,Y(J)/FLHT0
!     *			,EVEU(I,J),EVEV(I,J),FAI(I,J)
 !    *            ,EVETC(I,J)	      
            ENDDO
	          EVEETA(I)=EVEETA(I)/NTEVE
                EVEBED(I)=EVEBED(I)/NTEVE
	      ENDDO

         WRITE(61,22)
     *	'ZONE T="',T,'",F=POINT,I=',JM1-1,',J=',IM1-1
22	   FORMAT(A8,X,F8.4,X,A12,I5,X,A4,I5)
          DO I=2,IM1
	DO J=2,JM1
	  WRITE(61,121)(X(I)-X(IETA(2)))/FLHT0,Y(J)/FLHT0,EVETC(I,J)
	ENDDO
	ENDDO
         WRITE(61,22)
     *	'ZONE T="',T,'",F=POINT,I=',JM1-1,',J=',IM1-1
          DO I=2,IM1
	DO J=2,JM1
	  WRITE(61,121)(X(I)-X(IETA(2)))/FLHT0,Y(J)/FLHT0,EVETC1(I,J)
	ENDDO
	ENDDO
         WRITE(61,22)
     *	'ZONE T="',T,'",F=POINT,I=',JM1-1,',J=',IM1-1
          DO I=2,IM1
	DO J=2,JM1
	  WRITE(61,121)(X(I)-X(IETA(2)))/FLHT0,Y(J)/FLHT0,EVETC2(I,J)
	ENDDO
	ENDDO

 !          CALL TWODANYTC(EVEU,EVEV,EVETC,EVETCN,TWAVE)

  
   !      WRITE(62,21)'ZONE T="',T,'",F=POINT'       
  !       DO I = 1, IM1
  !        WRITE(62,112)(XU(I)-XU(IETM)),EVEETA(I),
  !   *	EVEBED(I),EVETC(I,11),EVETC(I,5),EVETC(I,15)
  !       ENDDO
112   FORMAT(1X,F15.3,2X,15(F10.5,2x))
21	FORMAT(A8,X,F8.4,X,A12)

	 
		 DO I=1,NNETA
		   DO J=2,JMAX
			SUMU(I,J)=SUMU(I,J)/NTEVE
			SUMVIS(I,J)=SUMVIS(I,J)/NTEVE
	        SUMTK(I,J)=SUMTK(I,J)/NTEVE
	        SUMTC(I,J)=SUMTC(I,J)/NTEVE
		   ENDDO
		  ENDDO

      QSUMS=0
      QSUME=0
	DO J=2,JMAX	 
	    QSUMS=QSUMS+EVEU(2,J)*DELTY(J)
	    QSUME=QSUME+EVEU(IM1-NSPE,J)*DELTY(J)
	    currentu1(j)=EVEU(IM1-NSPE,j)
	ENDDO

       DELQ1=STARTQIN-QSUMS

	do j=2,jmax	
	  if((y(j)-0.5*DELTY(J)).le.(FLHT0+EVEETA(2)))then 	      
         currentu(j)=CURRENTU(j)+DELQ1/(FLHT0+EVEETA(2))
	  ELSE
	   CURRENTU(J)=0.
	  ENDIF
	enddo

      DELQ2=STARTQIN-QSUME

	do j=2,jmax	
	  if((y(j)-0.5*DELTY(J)).le.(flht0))then		    
         currentu1(j)=currentu1(j)+DELQ2/(FLHT0+EVEETA(IETA(5)))
        ELSE
	   CURRENTU1(J)=0.
	  ENDIF
	enddo

	WRITE(112,121)T,QSUMS,QSUME

         WRITE(51,21)'ZONE T="',T,'",F=POINT'
         WRITE(52,21)'ZONE T="',T,'",F=POINT'
         WRITE(53,21)'ZONE T="',T,'",F=POINT'
         YC=FLHT0+EVEETA(IETA(3))
	   do j=2,jmax
	    write(51,121)y(j),currentu(j),currentu1(j),
     *		(sumu(i,j),i=1,nneta)
     	    write(52,121)y(j),(sumTK(i,j),i=1,nneta)
	    COMVIS=YC*EVEBED(IETA(3))*CNK*(1-Y(J)/YC)*(Y(J)/YC)
          write(53,121)y(j),COMVIS,(sumVIS(i,j),i=1,nneta)
	  ENDDO

          DO I=1,IM1
	     DO J=1,JM1
	       EVETCN(I,J)=EVETC(I,J)
	       EVETCN1(I,J)=EVETC1(I,J)
	       EVETCN2(I,J)=EVETC2(I,J)		   		   	      
		   EVEU(I,J)=0.
	       EVEV(I,J)=0.
	       EVETC(I,J)=0.
	       EVETC1(I,J)=0.
	       EVETC2(I,J)=0.
	     ENDDO
	       EVEETA(I)=0.
	       EVEBED(I)=0.
	    ENDDO

	     DO I=1,NNETA
		  DO J=2,JMAX
			SUMU(I,J)=0.
			SUMVIS(I,J)=0.
	        SUMTK(I,J)=0.
	        SUMTC(I,J)=0.
		   ENDDO
		  ENDDO
	    ENDIF   
	  ENDDO
       ENDIF
	ENDIF
	 	
 121  FORMAT(2X,F8.3,X,35(F12.8,x))
 113  format(x,8(f15.6,x))	 	
	end

	SUBROUTINE      anytc
	USE VOF2D
      IMPLICIT        REAL*8(A-H,O-Z)
!	INCLUDE         'vof.inc'
	REAL*8          TCEVE(MI),TEMPTC(MI,MJ)


	TCQN=TCQ
	TCEN=TCE
	TCE=0
	TCQ=0
	TCA1=0
	DO I=2,IM1
	TCEVE(I)=0
	DO J=2,jm1
	IF(F(I,J).GE.EMF)THEN
	TCEVE(I)=TCEVE(I)+TC(I,J)*DELTY(J)*DELTX(I)*F(I,J)
	TCA1=TCA1+TC(I,J)*DELTY(J)*DELTX(I)*F(I,J)
	ENDIF
	ENDDO
	ENDDO


 21   FORMAT(A8,X,F12.4,X,A12)
	DO I=2,IM1
	TCEVE(I)=TCEVE(I)/TCA1
	ENDDO

	IF(MOD(ICYCLE,20).EQ.0)THEN
	WRITE(115,21)'ZONE T="',T,'",F=POINT'
	DO I=2,IM1
	WRITE(115,101)X(I),TCEVE(I) 
	ENDDO
	ENDIF


	TCE1=0
	DO I=2,IM1
	TCE1=TCE1+X(I)*TCEVE(I)
	ENDDO
	TCE=TCE1

c-------------将中心点始终放在ieta（2）附近    
	if(itc.eq.0)THEN
	itc=2
	ELSE
      tcevemax=tceve(IETA(2))
	itceve=ieta(2)
	do i=IETA(2)-2,iETA(2)+5
	if(tceve(i).ge.tcevemax)then
	itceve=i
	tcevemax=tceve(i)
	endif
	enddo

	IF(ITCEVE.GE.IETA(2)+4)THEN
	Itc=0
	DO I=2,IM1
	DO J=2,JM1
	TEMPTC(I,J)=TC(I,J)
	TC(I,J)=0.
	ENDDO
	ENDDO

	DO I=2,IM1
	DO J=2,JM1
	II=I+(itceve-ieta(2))
	IF(II.LE.IM1.and.ii.ge.2)TC(I,J)=TEMPTC(II,J)
	ENDDO
	ENDDO 

	TCA1=0
	DO I=2,IM1
	TCEVE(I)=0
	DO J=2,jm1
	IF(F(I,J).GE.EMF)THEN
	TCEVE(I)=TCEVE(I)+TC(I,J)*DELTY(J)*DELTX(I)*F(I,J)
	TCA1=TCA1+TC(I,J)*DELTY(J)*DELTX(I)*F(I,J)
	ENDIF
	ENDDO
	ENDDO

	DO I=2,IM1
	TCEVE(I)=TCEVE(I)/TCA1
	ENDDO

	TCE=0
	DO I=2,IM1
	TCE=TCE+X(I)*TCEVE(I)
	ENDDO
	else
	itc=1
	endif

	ENDIF
c-------------

      DISE=(TCE1-TCEN)/(DELT)
      
	DO I=2,IM1
	TCQ=TCQ+(X(I)-TCE)**2*TCEVE(I)
	ENDDO

	DISP=0.5*(TCQ-TCQN)/(DELT)

	IF(Itc.NE.2)THEN
	WRITE(114,101)T,tca1,TCE,DISE,TCQ,DISP
	EVEDISP=EVEDISP+DISP
	IWRITE1=IWRITE1+1
	ENDIF

      DO N1=1,NSUM
      IF(T.GT.(TWAVE*N1-0.5*DELT).
     *		 AND.T.LT.(TWAVE*N1+0.5*DELT))THEN
	EVEDISP=EVEDISP/IWRITE1
	WRITE(116,101)T,EVEDISP
	EVEDISP=0
	iwrite1=0
	ENDIF
      ENDDO


101   FORMAT(2X,6(F12.6,3X))

      IF(TCEVE(IM1-nspe).GE.1.0D-4)THEN     
	WRITE(*,*)'X IS TOO SHORT'
	PAUSE
	ENDIF 

	END 

	SUBROUTINE      CALTRAGICY(FAI,UU,VV)
	USE VOF2D
      IMPLICIT        REAL*8(A-H,O-Z)
	REAL*8          FAI(MI,MJ),UU(MI,MJ),VV(MI,MJ)

	   DO I=1,IM1
	   DO J=1,JM1
	    FAI(I,J)=0.
	    M=1
	    N=1
	   DO  WHILE(M.LT.I.OR.N.LT.J)
	   IF(M.LT.I)THEN  
	    FAI(I,J)=FAI(I,J)-VV(M,N)*(X(M+1)-X(M))
	    M=M+1
	   ENDIF
	   IF(N.LT.J)THEN
	    FAI(I,J)=FAI(I,J)+UU(M,N)*(Y(N+1)-Y(N))
	    N=N+1
	   ENDIF
	   ENDDO

	  ENDDO
	  ENDDO

	END

	SUBROUTINE      TWODANYTC(UU,VV,TCC,TCCN,DT)
	USE             NUMERICAL_LIBRARIES
	USE VOF2D
      IMPLICIT        REAL*8(A-H,O-Z)
      PARAMETER       (NCP=41 )
      REAL*8          YANYTC(NCP),XANYTC(NCP,2),ANYTCB(2)
      REAL*8          UU(MI,MJ),VV(MI,MJ),TCC(MI,MJ),TCCN(MI,MJ)
         
	WRITE(301,21)'ZONE T="',T,'",F=POINT'
 21   FORMAT(A8,X,F12.4,X,A12) 
 
      IIS=IETA(2)
	IIE=IETA(2)+NCP-1

	DO J=2,JM1
       IF(YV(J-1).LT.FLHT0.AND.YV(J).GT.FLHT0)JJ=J
	ENDDO

      DO  J=2,JJ
	DO I=IIS,IIE
	    UEVE=0.5*(UU(I,J)+UU(I-1,J))
          VEVE=0.5*(VV(I,J)+VV(I,J-1))
          ADV=ADVECTION(3,I,J,EVEU,EVEV,TCC)

          RATIO1 = 0.5*DELTX(I+1)/DELTXU(I)
	    RATIO2 = 0.5*DELTX(I)/DELTXU(I-1)
          DDXE = (TCC(I+1,J)-TCC(I,J))/DELTXU(I)
          DDXW = (TCC(I,J)-TCC(I-1,J))/DELTXU(I-1)
          DXDX = (DDXE-DDXW)/DELTX(I)
	
C  -----CALCULATE DYDY
          RATIO1 = 0.5*DELTY(J+1)/DELTYV(J)
          RATIO2 = 0.5*DELTY(J)/DELTYV(J-1)
          DDYT= (TCC(I,J+1)-TCC(I,J))/DELTYV(J)
          DDYB= (TCC(I,J)-TCC(I,J-1))/DELTYV(J-1)
          DYDY= (DDYT-DDYB)/DELTY(J)
	    
		DEVETCDT=(TCC(I,J)-TCCN(I,J))/DT
	   	YANYTC(I-IIS+1)=DEVETCDT+ADV 

	    XANYTC(I-IIS+1,1)=DXDX
	    XANYTC(I-IIS+1,2)=DYDY
	  ENDDO
	    CALL DRLSE(NCP,YANYTC,2,XANYTC,NCP,0,ANYTCB,SST,SSE)
	    RANYTC=DSQRT(SST/(SST+SSE))

	    WRITE(301,101)Y(J)/FLHT0,ANYTCB(1),ANYTCB(2),SST,SSE,RANYTC
	ENDDO
101   FORMAT(X,8(F15.10,X))
    	
	END	  
	   
	SUBROUTINE  setpoint
	USE VOF2D
      real*8    x0,y0,tp,w0,v0,Dts,TTCC(3,MI,MJ),eyp1,eyp2

	dts=0.0064  !扩散系数
!	m0=0.01*uin*flht0
      ii=ieta(2)
      x0=x(ii)
	Y0=Y(JY(2))
	tp=abs(x(ii+1)-x0)/uin
	j=JY(2)			
	DO N=-4,4
      yp1=-(y(j)-y0+2.0*FLHT0*N)**2/(4.0*dts*tp)
	yp2=-(y(j)+y0+2.0*FLHT0*N)**2/(4.0*dts*tp)
	eyp1=exp(yp1)
	eyp2=exp(yp2)
	CN1=CN1+(eyp1+eyp2)   
	 ENDDO
	mp0=0.1 !uin*dsqrt(4.0*pi*tp*dts)/cn1 ! 
	 		    
      DO M=1,3
	Y0=Y(JY(M))  
	do i=2,im1
	tp=(x(i)-x0)/uin
	do j=2,jm1
	if(f(i,j).ge.emf.and.tp.gt.0)then !
	CN1=0.
	DO N=-4,4
	CN1=CN1+   
     *  ( Dexp(-(y(j)-y0+2.0*FLHT0*N)**2/(4.0*dts*tp))+
     *    Dexp(-(y(j)+y0+2.0*FLHT0*N)**2/(4.0*dts*tp))
     *  )
	 ENDDO
	CN2=(mp0/(uin*dsqrt(4.0*pi*tp*dts)))
	CN=CN1*CN2
	endif
	TTCC(M,I,J)=CN
	enddo
	enddo
	ENDDO

      WRITE(77,22)
     *	'ZONE T="',TP,'",F=POINT,I=',jm1-1,',J=',im1-1      
	DO I=2,IM1
	DO J=2,JM1
	write(77,10)(x(i)-X0)/flht0,y(j)/flht0,TTCC(1,I,J)
	ENDDO
	ENDDO

      WRITE(77,22)
     *	'ZONE T="',TP,'",F=POINT,I=',jm1-1,',J=',im1-1      
	DO I=2,IM1
	DO J=2,JM1
	write(77,10)(x(i)-X0)/flht0,y(j)/flht0,TTCC(2,I,J)
	ENDDO
	ENDDO

      WRITE(77,22)
     *	'ZONE T="',TP,'",F=POINT,I=',jm1-1,',J=',im1-1      
	DO I=2,IM1
	DO J=2,JM1
	write(77,10)(x(i)-X0)/flht0,y(j)/flht0,TTCC(3,I,J)
	ENDDO
	ENDDO

22	FORMAT(A8,X,F12.4,X,A12,I5,X,A4,I5)
10    format(x,6(f10.5,x))

      end

	SUBROUTINE  setpoints
	USE VOF2D
	real*8    dts(mj)

	WT=5.0
	WH=1.0
	ustar=0.083
	HBHS=0.85

	DTCM=FLHT0*USTAR*0.41
	DTWM=0.4*FLHT0*WH/WT
	DTCWM=DTCM+HBHS*DTWM


	DO  J = 2, jm1
	  yeta  =y(j)/flht0
	  dts(j) = DTCWM*yeta*(1.0-yeta)
	  if(dts(j).le.0)dts(j)=viscos
      ENDDO

	call subsetpoints(im1,jm1,x,y,uin,flht0,ieta(2),jy,dts,mp0)

	end

	SUBROUTINE  subsetpoints(mi,mj,x,y,uin,flht0,ii,jy,dts,mp0)
	real*8    x(mi),y(mj),dts(mj)
      real*8    x0,y0,tp,w0,v0,TTCC(3,MI,MJ),eyp1,eyp2,mp0,uin,flht0
	integer   jy(3)
      
	pi=3.1415926
!	dts=0.0064  !扩散系数
!	m0=0.01*uin*flht0

      x0=x(ii)
	Y0=Y(JY(2))
	tp=abs(x(ii+1)-x0)/uin
	j=JY(2)			
!	DO N=-4,4
  !    yp1=-(y(j)-y0+2.0*FLHT0*N)**2/(4.0*dts(j)*tp)
!	yp2=-(y(j)+y0+2.0*FLHT0*N)**2/(4.0*dts(j)*tp)
!	eyp1=exp(yp1)
!	eyp2=exp(yp2)
!	CN1=CN1+(eyp1+eyp2)   
!	 ENDDO
!	mp0=0.1 !uin*dsqrt(4.0*pi*tp*dts(j))/cn1 ! 
	 		    
      DO M=1,3
	Y0=Y(JY(M))  
	do i=2,mi
	tp=(x(i)-x0)/uin
	do j=2,mj
	if(tp.gt.0)then !
	CN1=0.
	DO N=-4,4
	CN1=CN1+   
     *  ( Dexp(-(y(j)-y0+2.0*FLHT0*N)**2/(4.0*dts(j)*tp))+
     *    Dexp(-(y(j)+y0+2.0*FLHT0*N)**2/(4.0*dts(j)*tp))
     *  )
	 ENDDO
	CN2=(mp0/(uin*dsqrt(4.0*pi*tp*dts(j))))
	CN=CN1*CN2
	endif
	TTCC(M,I,J)=CN
	enddo
	enddo
	ENDDO
	  
	  OPEN(77,FILE='2d points.DAT',STATUS='UNKNOWN')
	  WRITE(77,*)'TITLE="points"'
	  WRITE(77,*)'VARIABLES="X","y","c"' 

      WRITE(77,22)
     *	'ZONE T="',TP,'",F=POINT,I=',mj-1,',J=',mi-1      
	DO I=2,mi
	DO J=2,mj
	write(77,10)(x(i)-X0)/flht0,y(j)/flht0,TTCC(1,I,J)
	ENDDO
	ENDDO

      WRITE(77,22)
     *	'ZONE T="',TP,'",F=POINT,I=',mj-1,',J=',mi-1      
	DO I=2,mi
	DO J=2,mj
	write(77,10)(x(i)-X0)/flht0,y(j)/flht0,TTCC(2,I,J)
	ENDDO
	ENDDO

      WRITE(77,22)
     *	'ZONE T="',TP,'",F=POINT,I=',mj-1,',J=',mi-1      
	DO I=2,mi
	DO J=2,mj
	write(77,10)(x(i)-X0)/flht0,y(j)/flht0,TTCC(3,I,J)
	ENDDO
	ENDDO
	close(77)

22	FORMAT(A8,X,F12.4,X,A12,I5,X,A4,I5)
10    format(x,6(f10.5,x))

      end