|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？注册
x
我用老师给的程序，他的耗散率特别大，这样残差一直据高不下，残差能达到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
谢谢诸位的光临，请求你们指导！ |
|
IP卡
狗仔卡
发表于 2004-4-11 09:38:21
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
