耗散率狂大，残差一直据高不下，什么原因？

happyguy

我用老师给的程序，他的耗散率特别大，这样残差一直据高不下，残差能达到1.0E+8，耗散率通常在1.0E5左右，在壁面附近更高，源程序如下：
     ENTRY MODED
C----OUTLET
     DO 61 J=2,NJM1
     AE(NIM1,J)=0.0  
  61 AW(2,J)=0.0
C------TOP WALL
     J=2
     DO 62 I=nxin,NI-NXIN
     YP=0.5*SYCV(I,J)
     TERM=CMU75/(CAPPA*YP)
     SU(I,J)=GREAT*TERM*ABS(TE(I,J))**1.5
     SP(I,J)=-GREAT
AN(I,NJM1)=0
62  CONTINUE
C-----SYMMETRY
     DO 63 I=2,NXIN-1
     ED(I,1)=ED(I,2)      
ED(NI-I,1)=ED(NI-I,2)
     AS(NI-I,2)=0.0
AS(I,2)=0.0
AN(I,NJM1)=0
AN(NI-I,NJM1)=0.0
ED(I,2)=(ED(I,2)+ED(NI-I,2))/2.
ED(NI-I,2)=(ED(I,2)+ED(NI-I,2))/2.
63  CONTINUE
     RETURN

     DO 100 I=2,NI
     DO 100 J=2,NJ
C***********GEOMETRIC QUANTITIES USED TO CALCULATE CONVECTION
C***********AND DIFFUSION FLUXES THROUGH WEST AND SOUTH FACES OF
C***********E-CONSERVATION CONTROL VOLUMES
     AREAW=AEW(I,J)
     AREAS=ANS(I,J)
C----DISTANCES USED IN QUICK(ONLY WHEN NQUICK=1)
     IF(NQUICK.NE.1) GOTO 105
     DELX1=0.5*SXCV(I-1,J)
     DELX2=0.5*SXCV(I,J)
     DELY1=0.5*SYCV(I,J-1)
     DELY2=0.5*SYCV(I,J)
     B3X=DELX1*DELX2*(DELX1+DELX2)
     B3Y=DELY1*DELY2*(DELY1+DELY2)
105 CONTINUE
C----SPECIFICATION OF THE CONTROL VOLUME FACES WHERE EITHER PLDS
C----OR QUICK IS TO BE APPLIED,OUTSIDE THIS REGION THE CENTRAL
C----DIFFERENCING SCHEME IS USED
     LX=I.GT.2.AND.I.LT.NI
     LY=J.GT.2.AND.J.LT.NJ
C*********** CALCULATE CONVECTION AT WEST AND SOUTH FACES
     DENW=DEN(I-1,J)+FX(I-1,J)*(DEN(I,J)-DEN(I-1,J))
     DENS=DEN(I,J-1)+FY(I,J-1)*(DEN(I,J)-DEN(I,J-1))
     CW=DENW*U(I-1,J)*AREAW
     CS=DENS*V(I,J-1)*AREAS
C*********** CALCULATE DIFFUSION AT WEST AND SOUTH FACES
     GAMW=VIS(I-1,J)+FX(I-1,J)*(VIS(I,J)-VIS(I-1,J))
     GAMS=VIS(I,J-1)+FY(I,J-1)*(VIS(I,J)-VIS(I,J-1))
     GAMW=(GAMW-VISCOS)/PRED+VISCOS
     GAMS=(GAMS-VISCOS)/PRED+VISCOS
     DW=(GAMW+CE*VISXX(I-1,J))*AREAW/(0.5*(SXCV(I-1,J)+SXCV(I,J)))
     DS=(GAMS+CE*VISYY(I,J-1))*AREAS/(0.5*(SYCV(I,J-1)+SYCV(I,J)))
C*********** SCHEMES
     NSUP=-1
     IF(CW.GT.0.0)NSUP=1
     UP=0.5*FLOAT(1+NSUP)
     UM=0.5*FLOAT(1-NSUP)
     NSVP=-1
     IF(CS.GT.0.0)NSVP=1
     VP=0.5*FLOAT(1+NSVP)
     VM=0.5*FLOAT(1-NSVP)
     IF(LX.AND.NQUICK.EQ.1)GOTO 111
C----PLDS APPROXIMATION OF CONVECTION AT WEST FACE OF (I,J)
C----CONTROL VOLUME WHICH IS ALSO EAST FACE OF (I-1,J) C. V.
     TEMPW=DW-ABS(CW)*0.1
     COEW=0.
     TEMPW=TEMPW/(DW+TINY)
     IF(TEMPW.GT.TINY)COEW=DW*TEMPW**5
C---APPLICATION OF CENTRAL DIFFERENCING FOR THE FACES WHERE QUICK
C---AND PLDS ARE NOT USED
     IF(.NOT.LX)COEW=DW-MAX(FX(I-1,J)*CW,-(1.-FX(I-1,J))*CW)
C----SET SOURCES AND COEFFICIENTS DUE TO THE CONVECTION SCHEME TO
C----ZERO IN REGIONS WHERE QUICK IS NOT USED
     SPPX=0.
     SUPX=0.
     AE(I,J)=0.
     GOTO 300
C----QUICK APPROXIMATION OF CONVECTION AT WEST FACE OF (I,J) C.V.
C----WHICH IS ALSO EAST FACE OF (I-1,J) C.V.
111 KXP=I-1-NSUP
     LXP=KXP+(1-NSUP)/2
     DELX=0.5*(SXCV(I,J)+SXCV(I+1,J))
     IF(NSUP.EQ.1) DELX=0.5*(SXCV(I-1,J)+SXCV(I-2,J))
     DELX3=UM*DELX2-UP*DELX1-DELX*NSUP
     B1X=DELX2*DELX3*(DELX2-DELX3)
     B2X=DELX1*DELX3*(DELX1+DELX3)
     BX=B1X-B2X+B3X
     TEMP=B3X/BX*CW
     COEW=DW+(UM*B1X+UP*B2X)/BX*CW
     CW=CW-TEMP
     TEMPP=TEMP*UP
     TEMPM=TEMP*UM
C----SOURCES AND COEFF. DUE TO APPLYING QUICK TO WEST FACE OF
C----(I,J) C.V. AND EAST FACE OF (I-1,J) C. V.
     SPPX=-TEMPP
     SUPX=TEMPP*ED(LXP,J)
     AW(I-1,J)=AW(I-1,J)-TEMPP
     AE(I,J)=TEMPM
     SP(I-1,J)=SP(I-1,J)+TEMPM
     SU(I-1,J)=SU(I-1,J)-TEMPM*ED(LXP,J)
300 CONTINUE
     IF(LY.AND.NQUICK.EQ.1)GOTO 222
C----PLDS APPROXIMATION OF CONVECTION AT SOUTH FACE OF (I,J)
C----CONTROL VOLUME WHICH IS ALSO NORTH FACE OF (I,J-1) C. V.
     TEMPS=DS-ABS(CS)*0.1
     COES=0.
     TEMPS=TEMPS/(DS+TINY)
     IF(TEMPS.GT.TINY)COES=DS*TEMPS**5
C---APPLICATION OF CENTRAL DIFFERENCING FOR THE FACES WHERE QUICK
C---AND PLDS ARE NOT USED
     IF(.NOT.LY)COES=DS-MAX(FY(I,J-1)*CS,-(1.-FY(I,J-1))*CS)
C----SET SOURCES AND COEFFICIENTS DUE TO THE CONVECTION SCHEME TO
C----ZERO IN REGIONS WHERE QUICK IS NOT USED
     SPPY=0.
     SUPY=0.
     AN(I,J)=0.
     GOTO 400
C----QUICK APPROXIMATION OF CONVECTION AT SOUTH FACE OF (I,J) C.V.
C----WHICH IS ALSO NORTH FACE OF (I,J-1) C.V.
222 KYP=J-1-NSVP
     LYP=KYP+(1-NSVP)/2
     DELY=0.5*(SYCV(I,J)+SYCV(I,J+1))
     IF(NSVP.EQ.1) DELY=0.5*(SYCV(I,J-1)+SYCV(I,J-2))
     DELY3=VM*DELY2-VP*DELY1-DELY*NSVP
     B1Y=DELY2*DELY3*(DELY2-DELY3)
     B2Y=DELY1*DELY3*(DELY1+DELY3)
     BY=B1Y-B2Y+B3Y
     TEMP=B3Y/BY*CS
     COES=DS+(VM*B1Y+VP*B2Y)/BY*CS
     CS=CS-TEMP
     TEMPP=TEMP*VP
     TEMPM=TEMP*VM
C----SOURCES AND COEFF. DUE TO APPLYING QUICK TO SOUTH FACE OF
C----(I,J) C.V. AND NORTH FACE OF (I,J-1) C. V.
     SPPY=-TEMPP
     SUPY=TEMPP*ED(I,LYP)
     AS(I,J-1)=AS(I,J-1)-TEMPP
     AN(I,J)=TEMPM
     SP(I,J-1)=SP(I,J-1)+TEMPM
     SU(I,J-1)=SU(I,J-1)-TEMPM*ED(I,LYP)
400 CONTINUE
C*********** FALSE SOURCES
     CP=MAX(0.0,(SMPW(J)+CW))
     SMPW(J)=-CW-CS
     SMPW(J-1)=SMPW(J-1)+CS
     SU(I-1,J)=SU(I-1,J)+CP*ED(I-1,J)
     SP(I-1,J)=SP(I-1,J)-CP
C*********** ASSEMBLE COEFFICIENTS
     AW(I,J)=COEW+CW*UP
     AE(I-1,J)=AW(I,J)-CW+AE(I-1,J)
     AS(I,J)=COES+CS*VP
     AN(I,J-1)=AS(I,J)-CS+AN(I,J-1)
C----SOURCES DUE TO APPROXIMATION OF WEST AND SOUTH FACE valueS
C----USING QUICK
     SU(I,J)=SUPX+SUPY
     SP(I,J)=SPPX+SPPY
     IF(I.EQ.NI.OR.J.EQ.NJ) GOTO 100
C----DENSITY AT THE CENTER OF THE C.V.
     DENP=DEN(I,J)
C----VOLUME OF THE C.V.
     VOLP=VOL(I,J)
C********EVALUATION OF THE SOURCES DUE TO RESIDUAL DIFFUSION OF E
C----IN THE ASM EQUATIONS ONLY
     SUD=0.0
C******EVALUATION OF THE SOURCES IN THE DISSIPATION RATE EQUATION
C----K-E BVM
     SUC=C1*GEN(I,J)*ED(I,J)/(TE(I,J)+TINY)*VOLP
     TVIS=VIS(I,J)-VISCOS+TINY
     SPC=-C2*CMU*DEN(I,J)**2*ABS(TE(I,J))/TVIS*VOLP
     SU(I,J)=SU(I,J)+MAX(SUD,0.)+SUC
     SP(I,J)=SP(I,J)+MIN(SUD,0.)/(ED(I,J)+TINY)+SPC
100 CONTINUE
C*********** PROBLEM MODIFICATION
     CALL MODED
C****** FINAL COEFFICIENT ASSEMBLY AND RESIDUAL SOURCE CALCULATION
     RESORE=0.
     FAC=1.-URFE
     DO 200 I=2,NIM1
     DO 200 J=2,NJM1
     AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)-SP(I,J)
     RESOR=AN(I,J)*ED(I,J+1)+AS(I,J)*ED(I,J-1)+AE(I,J)*ED(I+1,J)
    +     +AW(I,J)*ED(I-1,J)-AP(I,J)*ED(I,J)+SU(I,J)
     IF(-SP(I,J).GE.GREAT)RESOR=0.
     RESORE=RESORE+ABS(RESOR)
C*********** UNDER-RELAXATION
     AP(I,J)=AP(I,J)/URFE
     SU(I,J)=SU(I,J)+FAC*AP(I,J)*ED(I,J)
200 CONTINUE
C********** SOLUTION OF DIFFERENCE EQUATIONS
     CALL LISOLV(2,2,NIM1,NJM1,ED,NSWPE)
     RETURN
     END