# 2D CFD code using SIMPLE algorithm

 Register Blogs Members List Search Today's Posts Mark Forums Read

 June 19, 2002, 16:13 2D CFD code using SIMPLE algorithm #1 bfan Guest   Posts: n/a Sponsored Links 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, 02:19 Re: 2D CFD code using SIMPLE algorithm #2 Free_and_SIMPLE Guest   Posts: n/a It's self explanatory!

 June 22, 2002, 22:01 Re: 2D CFD code using SIMPLE algorithm #4 bfan Guest   Posts: n/a my Email: binfan_y@yahoo.com do u have readme file about it?

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post CFD Main CFD Forum 17 January 3, 2017 18:09 ThomasTC Main CFD Forum 7 September 26, 2016 05:35 HectorRedal Main CFD Forum 8 June 13, 2011 22:03 Tareq Al-shaalan Main CFD Forum 10 June 12, 1999 23:27 Lans Main CFD Forum 13 October 27, 1998 11:20