|
[Sponsors] |
June 19, 2002, 17:13 |
2D CFD code using SIMPLE algorithm
|
#1 |
Guest
Posts: n/a
|
I know many people can have this code. Do you have document to explain it? If yes, pls give a copy.
LOGICAL LSTOP COMMON/CNTL/LSTOP *----------------------------------------------------------------------- CALL USER(1) CALL SETUP(1) CALL USER(2) 100 CALL USER(3) CALL USER(4) CALL USER(5) IF(LSTOP) STOP CALL SETUP(2) GOTO 100 END *================================================= ==================== SUBROUTINE DIFLOW *----------------------------------------------------------------------- COMMON/COEF/FLOW,DIFF,ACOF *----------------------------------------------------------------------- ACOF=DIFF IF(FLOW.EQ.0.) RETURN TEMP=DIFF-ABS(FLOW)*0.1 ACOF=0. IF(TEMP.LE.0.) RETURN TEMP=TEMP/DIFF ACOF=DIFF*TEMP**5 RETURN END *================================================= ==================== SUBROUTINE SOLVE *----------------------------------------------------------------------- COMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22, 22), + AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22, 22), + X(22),XU(22),XDIF(22),XCV(22),XCVS(22), + Y(22),YV(22),YDIF(22),YCV(22),YCVS(22), + YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22), + R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22), + DU(22,22),DV(22,22),FV(22),FVP(22), + FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22) *----------- Arrays U, V and PC may be absent in "SOLVE" ------------- DIMENSION U(22,22),V(22,22),PC(22,22) EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F( 1,1,3),PC(1,1)) *----------------------------------------------------------------------- CHARACTER*10 TITLE COMMON/A/TITLE(13) LOGICAL LSOLVE,LPRINT,LBLK ! &,LSTOP COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, + IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, + IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NT IMES(10),RHOCON *----------------------------------------------------------------------- ISTF=IST-1 JSTF=JST-1 IT1=L2+IST IT2=L3+IST JT1=M2+JST JT2=M3+JST *----------------------------------------------------------------------- DO 100 NT=1,NTIMES(NF) DO 100 N=NF,NF *----------------------------------------------------------------------- IF(.NOT.LBLK(NF)) GOTO 110 PT(ISTF)=0. QT(ISTF)=0. DO 120 I=IST,L2 BL=0. BLP=0. BLM=0. BLC=0. DO 130 J=JST,M2 BL=BL+AP(I,J) IF(J.NE.M2) BL=BL-AJP(I,J) IF(J.NE.JST) BL=BL-AJM(I,J) BLP=BLP+AIP(I,J) BLM=BLM+AIM(I,J) BLC=BLC+CON(I,J) & +AIP(I,J)*F(I+1,J,N) & +AIM(I,J)*F(I-1,J,N) & +AJP(I,J)*F(I,J+1,N) & +AJM(I,J)*F(I,J-1,N) & -AP(I,J)*F(I,J,N) 130 CONTINUE DENOM=BL-PT(I-1)*BLM IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E35 PT(I)=BLP/DENOM QT(I)=(BLC+BLM*QT(I-1))/DENOM 120 CONTINUE BL=0. DO 140 II=IST,L2 I=IT1-II BL=BL*PT(I)+QT(I) DO 140 J=JST,M2 140 F(I,J,N)=F(I,J,N)+BL *----------------------------------------------------------------------- PT(JSTF)=0. QT(JSTF)=0. DO 150 J=JST,M2 BL=0. BLP=0. BLM=0. BLC=0. DO 160 I=IST,L2 BL=BL+AP(I,J) IF(I.NE.L2) BL=BL-AIP(I,J) IF(I.NE.IST) BL=BL-AIM(I,J) BLP=BLP+AJP(I,J) BLM=BLM+AJM(I,J) BLC=BLC+CON(I,J) & +AIP(I,J)*F(I+1,J,N) & +AIM(I,J)*F(I-1,J,N) & +AJP(I,J)*F(I,J+1,N) & +AJM(I,J)*F(I,J-1,N) & -AP(I,J)*F(I,J,N) 160 CONTINUE DENOM=BL-PT(J-1)*BLM IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E+35 PT(J)=BLP/DENOM QT(J)=(BLC+BLM*QT(J-1))/DENOM 150 CONTINUE BL=0. DO 170 JJ=JST,M2 J=JT1-JJ BL=BL*PT(J)+QT(J) DO 170 I=IST,L2 170 F(I,J,N)=F(I,J,N)+BL 110 CONTINUE *----------------------------------------------------------------------- DO 180 J=JST,M2 PT(ISTF)=0. QT(ISTF)=F(ISTF,J,N) DO 190 I=IST,L2 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) PT(I)=AIP(I,J)/DENOM TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N) QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM 190 CONTINUE DO 200 II=IST,L2 I=IT1-II 200 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) 180 CONTINUE *----------------------------------------------------------------------- DO 210 JJ=JST,M3 J=JT2-JJ PT(ISTF)=0. QT(ISTF)=F(ISTF,J,N) DO 220 I=IST,L2 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) PT(I)=AIP(I,J)/DENOM TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N) QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM 220 CONTINUE DO 230 II=IST,L2 I=IT1-II 230 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) 210 CONTINUE *----------------------------------------------------------------------- DO 240 I=IST,L2 PT(JSTF)=0. QT(JSTF)=F(I,JSTF,N) DO 250 J=JST,M2 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) PT(J)=AJP(I,J)/DENOM TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM 250 CONTINUE DO 260 JJ=JST,M2 J=JT1-JJ 260 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) 240 CONTINUE *----------------------------------------------------------------------- DO 270 II=IST,L3 I=IT2-II PT(JSTF)=0. QT(JSTF)=F(I,JSTF,N) DO 280 J=JST,M2 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) PT(J)=AJP(I,J)/DENOM TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM 280 CONTINUE DO 290 JJ=JST,M2 J=JT1-JJ 290 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) 270 CONTINUE *----------------------------------------------------------------------- 100 CONTINUE DO 300 J=2,M2 DO 300 I=2,L2 CON(I,J)=0. AP(I,J)=0. 300 CONTINUE RETURN END *================================================= ==================== SUBROUTINE SETUP(K) *----------------------------------------------------------------------- COMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22, 22), + AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22, 22), + X(22),XU(22),XDIF(22),XCV(22),XCVS(22), + Y(22),YV(22),YDIF(22),YCV(22),YCVS(22), + YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22), + R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22), + DU(22,22),DV(22,22),FV(22),FVP(22), + FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22) *----------- Arrays U, V and PC may be absent in "SOLVE" ------------- DIMENSION U(22,22),V(22,22),PC(22,22) EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F( 1,1,3),PC(1,1)) *----------------------------------------------------------------------- CHARACTER*10 TITLE COMMON/A/TITLE(13) LOGICAL LSOLVE,LPRINT,LBLK,LSTOP COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, + IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, + IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NT IMES(10),RHOCON COMMON/CNTL/LSTOP COMMON/SORC/SMAX,SSUM COMMON/COEF/FLOW,DIFF,ACOF *----------------------------------------------------------------------- 10 FORMAT(/15X,' COMPUTATION IN CARTESIAN COORDINATES') 20 FORMAT(/15X,'COMPUTATION FOR AXISYMMETRIC SITUATION') 30 FORMAT(/15X,' COMPUTATION IN POLAR COORDINATES') 40 FORMAT(14X,40(1H*),/) *----------------------------------------------------------------------- GOTO (1000,2000) K *----------------------------------------------------------------------- * ENTRY SETUP1 1000 NFMAX=10 NP=11 NRHO=12 NGAM=13 LSTOP=.FALSE. DO 1010 I=1,10 LSOLVE(I)=.FALSE. NTIMES(I)=1 1010 LBLK(I)=.TRUE. DO 1020 I=1,13 LPRINT(I)=.FALSE. 1020 RELAX(I)=1. LAST=5 TIME=0. ITER=0 DT=1.E+10 IPREF=1 JPREF=1 RHOCON=1. *----------------------------------------------------------------------- L2=L1-1 L3=L2-1 M2=M1-1 M3=M2-1 X(1)=XU(2) DO 1030 I=2,L2 1030 X(I)=0.5*(XU(I)+XU(I+1)) X(L1)=XU(L1) Y(1)=YV(2) DO 1040 J=2,M2 1040 Y(J)=0.5*(YV(J+1)+YV(J)) Y(M1)=YV(M1) DO 1050 I=2,L1 1050 XDIF(I)=X(I)-X(I-1) DO 1060 I=2,L2 1060 XCV(I)=XU(I+1)-XU(I) DO 1070 I=3,L2 1070 XCVS(I)=XDIF(I) XCVS(L2)=XCVS(L2)+XDIF(L1) XCVS(3)=XCVS(3)+XDIF(2) DO 1080 I=3,L3 XCVI(I)=0.5*XCV(I) 1080 XCVIP(I)=XCVI(I) XCVIP(2)=XCV(2) XCVI(L2)=XCV(L2) DO 1090 J=2,M1 1090 YDIF(J)=Y(J)-Y(J-1) DO 1100 J=2,M2 1100 YCV(J)=YV(J+1)-YV(J) DO 1110 J=3,M2 1110 YCVS(J)=YDIF(J) YCVS(3)=YCVS(3)+YDIF(2) YCVS(M2)=YCVS(M2)+YDIF(M1) IF(MODE.NE.1) GOTO 1120 DO 1130 J=1,M1 RMN(J)=1.0 1130 R(J)=1.0 GOTO 1140 1120 DO 1150 J=2,M1 1150 R(J)=R(J-1)+YDIF(J) RMN(2)=R(1) DO 1160 J=3,M2 1160 RMN(J)=RMN(J-1)+YCV(J-1) RMN(M1)=R(M1) 1140 CONTINUE DO 1170 J=1,M1 SX(J)=1.0 SXMN(J)=1.0 IF(MODE.NE.3) GOTO 1170 SX(J)=R(J) IF(J.NE.1) SXMN(J)=RMN(J) 1170 CONTINUE DO 1180 J=2,M2 YCVR(J)=R(J)*YCV(J) ARX(J)=YCVR(J) IF(MODE.NE.3) GOTO 1180 ARX(J)=YCV(J) 1180 CONTINUE DO 1190 J=4,M3 1190 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J) YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3) YCVRS(M2)=.5*(R(M1)+R(M3))*YCVS(M2) IF(MODE.NE.2) GOTO 1200 DO 1210 J=3,M3 ARXJ(J)=.25*(1.+RMN(J)/R(J))*ARX(J) 1210 ARXJP(J)=ARX(J)-ARXJ(J) GOTO 1220 1200 DO 1230 J=3,M3 ARXJ(J)=.5*ARX(J) 1230 ARXJP(J)=ARXJ(J) 1220 ARXJP(2)=ARX(2) ARXJ(M2)=ARX(M2) DO 1240 J=3,M3 FV(J)=ARXJP(J)/ARX(J) 1240 FVP(J)=1.0-FV(J) DO 1250 I=3,L2 FX(I)=.5*XCV(I-1)/XDIF(I) 1250 FXM(I)=1.-FX(I) FX(2)=0.0 FXM(2)=1. FX(L1)=1. FXM(L1)=0. DO 1260 J=3,M2 FY(J)=.5*YCV(J-1)/YDIF(J) 1260 FYM(J)=1.-FY(J) FY(2)=0. FYM(2)=1.0 FY(M1)=1.0 FYM(M1)=0. *---------- CON,AP,U,V,RHO,PC,P ARRAYS ARE INITIALIZED HERE ---------- DO 1270 J=1,M1 DO 1270 I=1,L1 PC(I,J)=0. U(I,J)=0. V(I,J)=0. CON(I,J)=0. AP(I,J)=0. RHO(I,J)=RHOCON P(I,J)=0. 1270 CONTINUE IF(MODE.EQ.1) WRITE(10,10) IF(MODE.EQ.2) WRITE(10,20) IF(MODE.EQ.3) WRITE(10,30) WRITE(10,40) RETURN *----------------------------------------------------------------------- * ENTRY SETUP2 *----------------- COEFFICIENTS FOR THE U EQUATION ----------------- 2000 NF=1 IF(.NOT.LSOLVE(NF)) GOTO 2110 IST=3 JST=2 CALL USER(6) REL=1.-RELAX(NF) DO 2120 I=3,L2 FL=XCVI(I)*V(I,2)*RHO(I,1) FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1) FLOW=R(1)*(FL+FLM) DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2) CALL DIFLOW 2120 AJM(I,2)=ACOF+AMAX1(0.,FLOW) DO 2130 J=2,M2 FLOW=ARX(J)*U(2,J)*RHO(1,J) DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J)) CALL DIFLOW AIM(3,J)=ACOF+AMAX1(0.,FLOW) DO 2130 I=3,L2 IF(I.EQ.L2) GOTO 2140 FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J)) FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J) ) FLOW=ARX(J)*0.5*(FL+FLP) DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J)) GOTO 2150 2140 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J)) 2150 CALL DIFLOW AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) AIP(I,J)=AIM(I+1,J)-FLOW IF(J.EQ.M2) GOTO 2160 FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*R HO(I,J)) FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)* + RHO(I-1,J)) GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+ + 1.0E-30)*XCVI(I) GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)* + GAM(I-1,J)+1.E-30)*XCVIP(I-1) DIFF=RMN(J+1)*2.*(GM+GMM) GOTO 2170 2160 FL=XCVI(I)*V(I,M1)*RHO(I,M1) FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1) DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1) 2170 FLOW=RMN(J+1)*(FL+FLM) CALL DIFLOW AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) AJP(I,J)=AJM(I,J+1)-FLOW VOL=YCVR(J)*XCVS(I) APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1))/(XCVS(I)*DT) AP0/DX/DY AP(I,J)=AP(I,J)-APT AP->Sp CON(I,J)=CON(I,J)+APT*U(I,J) CON->Sc must be given firstly AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJM(I,J)+AJP(I,J))/ + RELAX(NF) CON(I,J)=CON(I,J)*VOL+AP(I,J)*U(I,J)*REL DU(I,J)=VOL/(XDIF(I)*SX(J)) CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J)) DU(I,J)=DU(I,J)/AP(I,J) 2130 CONTINUE CALL SOLVE 2110 CONTINUE *----------------- COEFFICIENTS FOR THE V EQUATION ----------------- NF=2 IF(.NOT.LSOLVE(NF)) GOTO 2210 IST=2 JST=3 CALL USER(6) REL=1.-RELAX(NF) DO 2220 I=2,L2 AREA=R(1)*XCV(I) FLOW=AREA*V(I,2)*RHO(I,1) DIFF=AREA*GAM(I,1)/YCV(2) CALL DIFLOW 2220 AJM(I,3)=ACOF+AMAX1(0.,FLOW) DO 2230 J=3,M2 FL=ARXJ(J)*U(2,J)*RHO(1,J) FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1) FLOW=FL+FLM DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J)) CALL DIFLOW AIM(2,J)=ACOF+AMAX1(0.,FLOW) DO 2230 I=2,L2 IF(I.EQ.L2) GOTO 2240 DIFF=2.*(GM+GMM)/SXMN(J) GOTO 2250 2240 FL=ARXJ(J)*U(L1,J)*RHO(L1,J) FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1) DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J)) 2250 FLOW=FL+FLM CALL DIFLOW AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) AIP(I,J)=AIM(I+1,J)-FLOW IF(J.EQ.M2) GOTO 2260 AREA=R(J)*XCV(I) FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J) FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J) )*RMN(J+1) FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I) DIFF=AREA*GAM(I,J)/YCV(J) GOTO 2270 2260 AREA=R(M1)*XCV(I) FLOW=AREA*V(I,M1)*RHO(I,M1) DIFF=AREA*GAM(I,M1)/YCV(M2) 2270 CALL DIFLOW AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) AJP(I,J)=AJM(I,J+1)-FLOW VOL=YCVRS(J)*XCV(I) SXT=SX(J) IF(J.EQ.M2) SXT=SX(M1) SXB=SX(J-1) IF(J.EQ.3) SXB=SX(1) APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)* + RHO(I,J-1)*0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT) AP(I,J)=AP(I,J)-APT CON(I,J)=CON(I,J)+APT*V(I,J) AP(I,J)=(-AP(I,J)*VOL+AIM(I,J)+AIP(I,J)+AJM(I,J)+AJP(I,J))/ + RELAX(NF) CON(I,J)=CON(I,J)*VOL+AP(I,J)*V(I,J)*REL DV(I,J)=VOL/YDIF(J) CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J)) DV(I,J)=DV(I,J)/AP(I,J) 2230 CONTINUE CALL SOLVE 2210 CONTINUE *--------- COEFFICIENT FOR THE PRESSURE CORRECTION EQUATION ---------- NF=3 IF(.NOT.LSOLVE(NF)) GOTO 2310 IST=2 JST=2 CALL USER(6) SMAX=0. SSUM=0. DO 2320 J=2,M2 DO 2320 I=2,L2 VOL=YCVR(J)*XCV(I) 2320 CON(I,J)=CON(I,J)*VOL DO 2330 I=2,L2 ARHO=R(1)*XCV(I)*RHO(I,1) CON(I,2)=CON(I,2)+ARHO*V(I,2) 2330 AJM(I,2)=0. DO 2340 J=2,M2 ARHO=ARX(J)*RHO(1,J) CON(2,J)=CON(2,J)+ARHO*U(2,J) AIM(2,J)=0. DO 2340 I=2,L2 IF(I.EQ.L2) GOTO 2350 ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) FLOW=ARHO*U(I+1,J) CON(I,J)=CON(I,J)-FLOW CON(I+1,J)=CON(I+1,J)+FLOW AIP(I,J)=ARHO*DU(I+1,J) AIM(I+1,J)=AIP(I,J) GOTO 2360 2350 ARHO=ARX(J)*RHO(L1,J) CON(I,J)=CON(I,J)-ARHO*U(L1,J) AIP(I,J)=0. 2360 IF(J.EQ.M2) GOTO 2370 ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)* RHO(I,J)) FLOW=ARHO*V(I,J+1) CON(I,J)=CON(I,J)-FLOW AJP(I,J)=ARHO*DV(I,J+1) AJM(I,J+1)=AJP(I,J) GOTO 2380 2370 ARHO=RMN(M1)*XCV(I)*RHO(I,M1) CON(I,J)=CON(I,J)-ARHO*V(I,M1) AJP(I,J)=0. 2380 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J) PC(I,J)=0. SMAX=AMAX1(SMAX,ABS(CON(I,J))) SSUM=SSUM+CON(I,J) 2340 CONTINUE CALL SOLVE *--------- COME HERE TO CORRECT THE PRESSURE AND VELOCITIES ---------- DO 2390 J=2,M2 DO 2390 I=2,L2 P(I,J)=P(I,J)+PC(I,J)*RELAX(NP) IF(I.NE.2)U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J)) IF(J.NE.2)V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J)) 2390 CONTINUE 2310 CONTINUE *----------------- COEFFICIENTS FOR OTHER EQUATIONS ------------------ IST=2 JST=2 DO 2410 N=4,NFMAX NF=N IF(.NOT.LSOLVE(NF)) GOTO 2410 CALL USER(6) REL=1.-RELAX(NF) DO 2420 I=2,L2 AREA=R(1)*XCV(I) FLOW=AREA*V(I,2)*RHO(I,1) DIFF=AREA*GAM(I,1)/YDIF(2) CALL DIFLOW 2420 AJM(I,2)=ACOF+AMAX1(0.,FLOW) DO 2430 J=2,M2 FLOW=ARX(J)*U(2,J)*RHO(1,J) DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J)) CALL DIFLOW AIM(2,J)=ACOF+AMAX1(0.,FLOW) DO 2430 I=2,L2 IF(I.EQ.L2) GOTO 2440 FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)* RHO(I,J)) DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+ + XCV(I+1)*GAM(I,J)+1.E-30)*SX(J)) GOTO 2450 2440 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J)) 2450 CALL DIFLOW AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) AIP(I,J)=AIM(I+1,J)-FLOW AREA=RMN(J+1)*XCV(I) IF(J.EQ.M2) GOTO 2460 FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RH O(I,J)) DIFF=AREA*2.*GAM(I,J+1)*GAM(I,J)/(YCV(J)*GAM(I,J+1)+ + YCV(J+1)*GAM(I,J)+1.E-30) GOTO 2470 2460 FLOW=AREA*V(I,M1)*RHO(I,M1) DIFF=AREA*GAM(I,M1)/YDIF(M1) 2470 CALL DIFLOW AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) AJP(I,J)=AJM(I,J+1)-FLOW VOL=YCVR(J)*XCV(I) APT=RHO(I,J)/DT AP(I,J)=AP(I,J)-APT CON(I,J)=CON(I,J)+APT*F(I,J,NF) AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))/ + RELAX(NF) CON(I,J)=CON(I,J)*VOL+AP(I,J)*F(I,J,NF)*REL 2430 CONTINUE CALL SOLVE 2410 CONTINUE TIME=TIME+DT ITER=ITER+1 IF(ITER.EQ.LAST) LSTOP=.TRUE. RETURN END *================================================= ==================== SUBROUTINE SUPPLY(K) *----------------------------------------------------------------------- COMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22, 22), + AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22, 22), + X(22),XU(22),XDIF(22),XCV(22),XCVS(22), + Y(22),YV(22),YDIF(22),YCV(22),YCVS(22), + YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22), + R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22), + DU(22,22),DV(22,22),FV(22),FVP(22), + FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22) *----------- Arrays U, V and PC may be absent in "SOLVE" ------------- DIMENSION U(22,22),V(22,22),PC(22,22) EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F( 1,1,3),PC(1,1)) *----------------------------------------------------------------------- CHARACTER*10 TITLE COMMON/A/TITLE(13) LOGICAL LSOLVE,LPRINT,LBLK ! &,LSTOP COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, + IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, + IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NT IMES(10),RHOCON *----------------------------------------------------------------------- 110 FORMAT(1X,26(1H*),3X,A10,3X,26(1H*)) 120 FORMAT(1X,'I =',I8,6I10) 130 FORMAT(1X,'J') 140 FORMAT(1X,I2,2X,1P7E10.2) 150 FORMAT(1X,' ') 160 FORMAT(1X,'I =',I8,6I10) 170 FORMAT(1X,'X = ',1P7E10.2) 180 FORMAT(1X,'XU =',1P7E10.2) 190 FORMAT(1X,'TH =',1P7E10.2) 200 FORMAT(1X,'J =',I8,6I10) 210 FORMAT(1X,'Y = ',1P7E10.2) 220 FORMAT(1X,'YV =',1P7E10.2) *----------------------------------------------------------------------- GOTO (1000,2000) K *----------------------------------------------------------------------- * ENTRY UGRID 1000 XU(2)=0. DX=XL/FLOAT(L1-2) DO 1010 I=3,L1 1010 XU(I)=XU(I-1)+DX YV(2)=0. DY=YL/FLOAT(M1-2) DO 1020 J=3,M1 1020 YV(J)=YV(J-1)+DY RETURN *----------------------------------------------------------------------- * ENTRY PRINT 2000 IF(.NOT.LPRINT(3)) GOTO 2010 *------------------ CALCULATE THE STREAM FUNCTION ------------------ F(2,2,3)=0. DO 2020 I=2,L1 IF(I.NE.2) F(I,2,3)=F(I-1,2,3)-RHO(I-1,1)*V(I-1,2)*R(1)*XCV(I-1) DO 2020 J=3,M1 RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1) 2020 F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1) 2010 CONTINUE IF(.NOT.LPRINT(NP)) GOTO 2030 *----------- CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION ----------- DO 2040 J=2,M2 P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3) 2040 P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2) DO 2050 I=2,L2 P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3) 2050 P(I,M1)=(P(I,M2)*YCVS(M2)-P(I,M3)*YDIF(M1))/YDIF(M2) P(1,1)=P(2,1)+P(1,2)-P(2,2) P(L1,1)=P(L2,1)+P(L1,2)-P(L2,2) P(1,M1)=P(2,M1)+P(1,M2)-P(2,M2) P(L1,M1)=P(L2,M1)+P(L1,M2)-P(L2,M2) PREF=P(IPREF,JPREF) DO 2060 J=1,M1 DO 2060 I=1,L1 2060 P(I,J)=P(I,J)-PREF 2030 CONTINUE *----------------------------------------------------------------------- WRITE(10,150) IEND=1 2080 IF(IEND.EQ.L1) GOTO 2070 IBEG=IEND+1 IEND=IEND+7 IEND=MIN0(IEND,L1) WRITE(10,160) (I,I=IBEG,IEND) WRITE(10,180) (XU(I),I=IBEG,IEND) GOTO 2080 2070 WRITE(10,150) JEND=1 2100 IF(JEND.EQ.M1) GOTO 2090 JBEG=JEND+1 JEND=JEND+7 JEND=MIN0(JEND,M1) WRITE(10,200) (J,J=JBEG,JEND) WRITE(10,220) (YV(J),J=JBEG,JEND) GOTO 2100 2090 CONTINUE *----------------------------------------------------------------------- WRITE(10,150) IEND=0 2140 IF(IEND.EQ.L1) GOTO 2110 IBEG=IEND+1 IEND=IEND+7 IEND=MIN0(IEND,L1) WRITE(10,160) (I,I=IBEG,IEND) IF(MODE.EQ.3) GOTO 2120 WRITE(10,170) (X(I),I=IBEG,IEND) GOTO 2130 2120 WRITE(10,190) (X(I),I=IBEG,IEND) 2130 GOTO 2140 2110 WRITE(10,150) JEND=0 2160 IF(JEND.EQ.M1) GOTO 2150 JBEG=JEND+1 JEND=JEND+7 JEND=MIN0(JEND,M1) WRITE(10,200) (J,J=JBEG,JEND) WRITE(10,210) (Y(J),J=JBEG,JEND) GOTO 2160 2150 CONTINUE *----------------------------------------------------------------------- DO 2170 N=1,NGAM NF=N IF(.NOT.LPRINT(NF)) GOTO 2170 WRITE(10,150) WRITE(10,110) TITLE(NF) IFST=1 JFST=1 IF(NF.EQ.1.OR.NF.EQ.3) IFST=2 IF(NF.EQ.2.OR.NF.EQ.3) JFST=2 IBEG=IFST-7 2190 CONTINUE IBEG=IBEG+7 IEND=IBEG+6 IEND=MIN0(IEND,L1) WRITE(10,150) WRITE(10,120) (I,I=IBEG,IEND) WRITE(10,130) JFL=JFST+M1 DO 2180 JJ=JFST,M1 J=JFL-JJ WRITE(10,140) J,(F(I,J,NF),I=IBEG,IEND) 2180 CONTINUE IF(IEND.LT.L1) GOTO 2190 2170 CONTINUE RETURN END |
|
June 20, 2002, 03:19 |
Re: 2D CFD code using SIMPLE algorithm
|
#2 |
Guest
Posts: n/a
|
It's self explanatory!
|
|
June 22, 2002, 12:50 |
Re: 2D CFD code using SIMPLE algorithm
|
#3 |
Guest
Posts: n/a
|
please contact me i try to help you.
|
|
June 22, 2002, 23:01 |
Re: 2D CFD code using SIMPLE algorithm
|
#4 |
Guest
Posts: n/a
|
my Email: binfan_y@yahoo.com do u have readme file about it?
|
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
CFD Salary | CFD | Main CFD Forum | 17 | January 3, 2017 18:09 |
Looking for simple, robust interface tracking algorithm | ThomasTC | Main CFD Forum | 7 | September 26, 2016 06:35 |
Has enyone developped code for CFD using FEM? | HectorRedal | Main CFD Forum | 8 | June 13, 2011 23:03 |
Which is better to develop in-house CFD code or to buy a available CFD package. | Tareq Al-shaalan | Main CFD Forum | 10 | June 13, 1999 00:27 |
What kind of Cmmercial CFD code you feel well? | Lans | Main CFD Forum | 13 | October 27, 1998 11:20 |