CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

2D CFD code using SIMPLE algorithm

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

Reply
 
LinkBack Thread Tools Display Modes
Old   June 19, 2002, 16:13
Default 2D CFD code using SIMPLE algorithm
  #1
bfan
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
  Reply With Quote

Old   June 20, 2002, 02:19
Default Re: 2D CFD code using SIMPLE algorithm
  #2
Free_and_SIMPLE
Guest
 
Posts: n/a
It's self explanatory!
  Reply With Quote

Old   June 22, 2002, 11:50
Default Re: 2D CFD code using SIMPLE algorithm
  #3
moosa
Guest
 
Posts: n/a
please contact me i try to help you.
  Reply With Quote

Old   June 22, 2002, 22:01
Default 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?
  Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Looking for simple, robust interface tracking algorithm ThomasTC Main CFD Forum 6 February 17, 2012 19:13
Has enyone developped code for CFD using FEM? HectorRedal Main CFD Forum 8 June 13, 2011 22:03
CFD Salary CFD Main CFD Forum 15 September 4, 1999 14:04
Which is better to develop in-house CFD code or to buy a available CFD package. Tareq Al-shaalan Main CFD Forum 10 June 12, 1999 23:27
What kind of Cmmercial CFD code you feel well? Lans Main CFD Forum 13 October 27, 1998 11:20


All times are GMT -4. The time now is 08:38.