C**************************************************************** C NATURAL CONVECTION IN NON-NEWTONIAN CAVITY C Younis et al case C**************************************************************** BLOCK DATA CHARACTER*10 TITLE LOGICAL LSOLVE,LPRINT,LBLK,LSTOP COMMON/CHAR/TITLE(13) COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, 1 IST,JST,ITER,LAST,RELAX(13),TTIME,DT,XL,YL, 2 IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON COMMON/CNTL/LSTOP COMMON/SORC/SMAX,SSUM COMMON/COEF/FLOW,DIFF,ACOF DATA NFMAX,NP,NRHO,NGAM/10,11,12,13/ DATA LSTOP/1*.FALSE./ DATA LBLK/10*.TRUE./ DATA MODE,TTIME,ITER/1,0.,0/ DATA RELAX(13)/1./ DATA NTIMES/10*1/ DATA DT,IPREF,JPREF,RHOCON/1.E+10,1,1,1./ C**************************************************************** C ENTRADA DE DATOS A ETIQUETA COMUN USANDO DATOS DECLARADOS C *********************************************************** DATA TITLE(1),TITLE(2),TITLE(4),TITLE(5),TITLE(11)/ 1'VEL U',' VEL V',' TEMP','COn','PRESSURE'/ DATA RELAX(1),RELAX(2),RELAX(4),RELAX(5),RELAX(11) 1/0.25,0.25,0.5,0.5,0.8/ DATA (LSOLVE(I),LPRINT(I),I=1,5)/10*.TRUE./ DATA LAST/300000000/ END C**************************************************************** C DEFINICIÓN DE LOS ARCHIVOS DE DATOS C**************************************************************** SUBROUTINE OPENF OPEN(UNIT=1,FILE='DATOS.TXT') OPEN(UNIT=1504,FILE='Todo_TpoFinal.PLT') OPEN(UNIT=1505,FILE='TodoAdmim.PLT') OPEN(UNIT=1506,FILE='TodoDimensional.PLT') OPEN(UNIT=73,FILE='TPROM.PLT') OPEN(UNIT=101,FILE='DURACION.TXT') RETURN END C**************************************************************** SUBROUTINE USER LOGICAL LSOLVE,LPRINT,LBLK,LSTOP PARAMETER (NOD=180) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), 1CON(NOD,NOD), 1 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), 2 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), 3 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), 4 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), 5 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD) COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, 1 IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, 2 IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCO COMMON/CNTL/LSTOP COMMON/SORC/SMAX,SSUM COMMON/COEF/FLOW,DIFF,ACOF DIMENSION U(NOD,NOD),V(NOD,NOD),ABC1(NOD,NOD),ABC2(NOD,NOD), 1EFG1(NOD,NOD),EFG2(NOD,NOD),CONVER1(NOD,NOD),CONVER2(NOD,NOD), 2GHI1(NOD,NOD),GHI2(NOD,NOD),CONVER3(NOD,NOD), 3JKL1(NOD,NOD),JKL2(NOD,NOD),CONVER4(NOD,NOD), 4UANT(NOD,NOD),VANT(NOD,NOD),TANT(NOD,NOD),CANT(NOD,NOD) EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),STR(1,1)) DIMENSION TH(NOD),THU(NOD),THDIF(NOD),THCV(NOD),THCVS(NOD) EQUIVALENCE(X,TH),(XU,THU),(XDIF,THDIF),(XCV,THCV), 1(XCVS,THCVS),(XL,THL) COMMON/MVAT/NTYV,NTXV,VTYV(30),VTXV(30),YLTV(30),XLTV(30),ITRAM, 1COORD,RX(30),RY(30) C**************************************************************** DIMENSION T(NOD,NOD),C(NOD,NOD),STR(NOD,NOD),RA(NOD,NOD), 1PR(NOD,NOD),dE(NOD,NOD) EQUIVALENCE(F(1,1,4),T(1,1)),(F(1,1,5),C(1,1)) DIMENSION ANULH(NOD),ANULC(NOD) DIMENSION ANULI(NOD),ANULD(NOD),ANULS(NOD),ANULB(NOD) DIMENSION ETA(NOD,NOD),ETAI(NOD,NOD),ETAE(NOD,NOD) REAL NUSLH,NUSLC,NUSMH,NUSMC,NUMINH,NUMAXH,RA0,PR0 real tdim(nod,nod),RHOvelu(nod,nod),RHOvelv(nod,nod) &,RHOtem(nod,nod),RHOconc(nod,nod),rhocp(nod,nod) REAL TA(2) REAL VCONTROL(NOD) integer PARTES(NOD,NOD) REAL VINIC,NVC,TFILM,tpromDim,visc,viscK,etain,varaux1 &,varaux2,varaux3 REAL xlong,ylong,tdimen,lo,xmax1,xmin2,ymax1,ymin2 INTEGER CONTI,TINIC,TEX,xnodo,ynodo,salida,imprimir1 C**************************************************************** C**************************************************************** ENTRY GRID MODE=1. L1 = 50 XL = 1 M1 = 50 YL = 1 call ugrid RETURN C**************************************************************** C**************************************************************** ENTRY START DTPO = 1.e-2 !time step TIEMPO = 0. TFIN = 1200 !TIEMPO FINAL TIEMPOFINAL=TIMEF() !TIEMPO COMPUTACIONAL EMAX=0. ETMAX=0. ITERMIN=7.;ITERMAX=150. !ITERACIÓN MÍNIMA Y MAXIMA CRITERIO=1.E-6 !criterio de convergencia para paso de tiempo TPRINT=0. TOUT=1e-2 !*************************************************************** salida = 0 imprime = 0 SMALL=1.E-30 !*****posición de los nodos DO J=1,M1 DO I=1,L1 IF(X(I).LE.0.05) xmax1 = i+1 If(x(i).le.1.-0.05) xmin2 = i-1 if(y(j).le.0.05) ymax1 = j if(y(j).le.1-0.05) ymin2 = j-1 ENDDO ENDDO partes = 0 DO J=1,M1 DO I=1,L1 IF(i.le.xmax1.or.i.gt.xmin2 & .or.j.le.ymax1.or.j.gt.ymin2)then partes(i,j) = 1 !es el sólido endif ENDDO ENDDO !***********************LIMITES DE TEMPERAURA FLUIDO**************************** TINIC = 45 !TEMPERATURA INICIAL GOLDEN SYRUP TEX = 25 !TEMPERATURA FINAL TFILM = (TEX+TINIC)/2. !Temperatura media !**************************lARGO CARACTERISTICO************************************* LO = 0.2 !LONGITUD dimensional R0 = 0.2 !Longitud dimensional GE = 9.81 !GRAVEDAD CTEGAS = 8.31451 !GASES J/mol K ALFAS = 1.0043E-7 !DIFUSIVIDAD TERM plexiglas !PARAMETROS MODELO DE VISCOSIDAD GOLDEN SYRUP******************** AA = -7.5907E-7 BB = 3.8968E-4 CC = 4.0130E-2 ETA0 = 4.485E-8 !VISCOSIDAD INICIAL PA S MODELO GOLDEN SY C**************** liquido interior ********************** DENS = 1.438E3 !DENSIDAD EN KG/M3 BETA = 4.33E-4 !EXP. VOLUMETRICA K-1 ETAin = ETA0*(EXP(1./((AA*TFILM**2)+ 1(BB*TFILM)+CC))) ETAREF=ETAin !VISCOCIDAD aparente inicial PA S U0=SQRT(GE*BETA*r0*(TINIC-TEX)) !MAXIMA VELOCIDAD EN EL FLUIDO ALFA = 1.21E-7 !DIFUSIVIDAD TERMICA GOLDEN SYRUP M2/S PR0 = ETAREF/dens/ALFA RA0=(BETA*GE*(TINIC-TEX)*R0**3) !NUMERO DE RAYLEIGH INICIAL &/((ETAREF/dens)*ALFA) !VALORES INICIALES y CONDICIONES DE BORDE INICIALES*********** DO 100 J=1,M1 DO 100 I=1,L1 U(I,J) = 0. V(I,J) = 0. t(i,j) = 0. T(I,1) = -1. !CB. ARRIBA T(I,M1) = -1. !CB. ARRIBA T(1,J) = -1. !CB. IZQUIERDA T(L1,J) = -1. !CB. DERECHA ETA(I,J) = 1. ETAI(I,J) = 1. ETAE(I,J) = 1. !INICIALIZACIÓN DE RA Y PR *********************************** RA(I,J) = RA0 PR(I,J) = PR0 100 CONTINUE !**************************************************************** DO 110 J=1,M1 DO 110 I=1,L1 UANT(I,J)=U(I,J) VANT(I,J)=V(I,J) TANT(I,J)=T(I,J) CANT(I,J)=C(I,J) ABC1(I,J)=1. EFG1(I,J)=1. GHI1(I,J)=1. JKL1(I,J)=1 110 CONTINUE RETURN C**************************************************************** C**************************************************************** ENTRY DENSE DO I=1,L1 DO J=1,M1 !solo para u y v velocidades if(partes(i,j).eq.1)rho(i,j) = 1. !solido if(partes(i,j).eq.0)then !fluido RHO(I,J) = SQRT(RA(I,J)/PR(I,J)) endif enddo enddo RETURN C**************************************************************** C**************************************************************** C**************************************************************** ENTRY BOUND DO 300 I=1,L1 DO 300 J=1,M1 T(I,1) = -1 !CB. ABAJO DEBE CORRESPONDER A 25 °C T(I,M1) = -1. !CB. ARRIBA T(1,J) = -1. !CB. IZQUIERDA T(L1,J) = -1. !CB. DERECHA !VELOCIDADES CERO EN LAS PAREDES IF(partes(I,J).EQ.1)U(I,J)=0 IF(partes(I,J).EQ.1)V(I,J)=0 300 CONTINUE RETURN C**************************************************************** C**************************************************************** C**************************************************************** ENTRY OUTPUT CONT=0 DO 521 I=1,L1 DO 521 J=1,M1 ABC2(I,J)=U(I,J) EFG2(I,J)=V(I,J) GHI2(I,J)=T(I,J) JKL2(I,J)=C(I,J) CONVER1(I,J)=ABS(ABC2(I,J)-ABC1(I,J)) CONVER2(I,J)=ABS(EFG2(I,J)-EFG1(I,J)) CONVER3(I,J)=ABS(GHI2(I,J)-GHI1(I,J)) CONVER4(I,J)=ABS(JKL2(I,J)-JKL1(I,J)) ABC1(I,J)=ABC2(I,J) EFG1(I,J)=EFG2(I,J) GHI1(I,J)=GHI2(I,J) JKL1(I,J)=JKL2(I,J) !criterio de convergencia para impresión de resultados y proximo paso de tpo IF(CONVER1(I,J).LE.CRITERIO.AND.CONVER2(I,J).LE.CRITERIO.AND. 1CONVER3(I,J).LE.CRITERIO.AND.CONVER4(I,J).LE.CRITERIO)THEN CONT=CONT+1 ENDIF 521 CONTINUE 333 FORMAT(1I8,F10.4,E14.6,F10.1,3E12.4) WRITE(*,*)ITERN,TIEMPO,tdimen,ETMAX,dtpo !INICIA IMPRESIÓN DE RESULTADOS Y PASO DE TIEMPO IF (CONT.GE.L1*M1*0.99.AND.ITERN.GE.ITERMIN.AND.SMAX.LE.1.E-4 !INTRODUCIR CRITERIO DE CONVERGENCIA %.OR.ITERN.EQ.itermax)THEN C*******variación maxima de la variable para aumentar el paso de tiempo en los primeros segundos ETMAX=0. DO 406 I=1,L1 DO 406 J=1,M1 ! IF(ABS(U(I,J)-UANT(I,J)).GT.ETMAX)ETMAX=ABS(U(I,J)-UANT(I,J)) IF(ABS(V(I,J)-VANT(I,J)).GT.ETMAX)ETMAX=ABS(V(I,J)-VANT(I,J)) ! IF(ABS(T(I,J)-TANT(I,J)).GT.ETMAX)ETMAX=ABS(T(I,J)-TANT(I,J)) ! IF(ABS(C(I,J)-CANT(I,J)).GT.ETMAX)ETMAX=ABS(C(I,J)-CANT(I,J)) 406 CONTINUE IF(ETMAX.LT.1E-4.AND.DTPO.LT.10)DTPO=DTPO*1.2 tdimen = (tiempo * LO)/U0 ! salidas a 89.2562, 1190.083 y 1785.12 segundos en modelo de younis varaux1 = 89.2562 varaux2 = 1190.083 varaux3 = 1785.12 IF(tdimen.GE.0.AND.salida.EQ.0) imprimir1 = 1 IF(tdimen.GE.VARAUX1.AND.salida.EQ.1) imprimir1 = 1 IF(tdimen.GE.VARAUX2.AND.salida.eq.2) imprimir1 = 1 IF(tdimen.GE.VARAUX3.AND.salida.eq.3) imprimir1 = 1 810 FORMAT (10F17.4) if(imprimir1.eq.1)then imprimir1 = 0 salida = salida + 1 WRITE(1505,*) 'TITLE = "TodoAdim"' WRITE(1505,*) 'VARIABLES="X","Y","U","V","T","C"' WRITE(1505,*) 'ZONE T="',TIEMPO,'" I=',L1,', J=',M1 DO 850 J=1,M1 DO 850 I=1,L1 WRITE(1505,810) X(I),Y(J),U(I,J),V(I,J),T(I,J),C(I,J) 850 CONTINUE WRITE(1506,*) 'TITLE = "TodoDim"' WRITE(1506,*) 'VARIABLES="X","Y","U","V","T","C"' WRITE(1506,*) 'ZONE T="',Tdimen,'" I=',L1,', J=',M1 DO 840 J=1,M1 DO 840 I=1,L1 WRITE(1506,810) X(I)*R0,Y(J)*R0,U0*U(I,J),U0*V(I,J), 1 ((TINIC-TEX)*T(I,J))+TINIC,C(I,J) 840 CONTINUE endif IF(TIEMPO.GE.TFIN+DTPO)THEN WRITE(1504,*) 'TITLE = "prueba1"' WRITE(1504,*) 'VARIABLES="X","Y","U","V","T","C"' WRITE(1504,*) 'ZONE T="',TIEMPO,'" I=',L1,', J=',M1 DO J=1,M1 DO I=1,L1 WRITE(1504,*) X(I),Y(J),U(I,J),V(I,J),T(I,J),C(I,J) enddo enddo DTIEMPO=ETIME(TA) WRITE(*,*) 'El programa funciono por:',DTIEMPO,'CPU segundos.' WRITE(*,*) 'DONDE:' WRITE(*,*) 'El tiempo de usuario fue',TA(1),'segundos' WRITE(*,*) 'Y el tiempo de sistema fue:',TA(2),'segundos' WRITE(*,*) 'El tiempo de usuario fue',TA(1)/60,'MINUTOS' WRITE(*,*) 'Y el tiempo de sistema fue:',TA(2)/60,'MINUTOS' WRITE(*,*) 'La cantidad total de iteraciones fue:',ITER WRITE(101,*) 'El programa funciono por:',DTIEMPO,'CPU segundos.' WRITE(101,*) 'DONDE:' WRITE(101,*) 'El tiempo de usuario fue',TA(1),'segundos' WRITE(101,*) 'Y el tiempo de sistema fue:',TA(2),'segundos' WRITE(101,*) 'El tiempo de usuario fue',TA(1)/60,'MINUTOS' WRITE(101,*) 'Y el tiempo de sistema fue:',TA(2)/60,'MINUTOS' PAUSE STOP ENDIF DO 470 I=1,L1 DO 470 J=1,M1 UANT(I,J)=U(I,J) VANT(I,J)=V(I,J) TANT(I,J)=T(I,J) cANT(I,J)=T(I,J) 470 CONTINUE !*****calculo de la temperatura promedio al interior de la cavidad conti = 0 TPROM = 0 DO 703 J=1,M1 DO 703 I=1,L1 if(partes(i,j).eq.0)then TPROM=T(I,J)+TPROM CONTI = CONTI + 1 !NUMERO DE VOLUMENES DE CONTROL CAMARA endif 703 CONTINUE tprom = tprom/conti WRITE(73,*) tdimen,tprom !la iteración vuelve a cero ITERN=0. !******aumento del paso de tiempo TIEMPO=TIEMPO+DTPO ENDIF ITERN=ITERN+1 RETURN C**************************************************************** C**************************************************************** ENTRY GAMSOR tdim = 0. etai = 0. ra = 0. pr = 0. DO 502 I=1,L1 DO 502 J=1,M1 tdim(i,j) = TINIC+T(I,J)*(TINIC-TEX) ETAI(I,J) = (ETA0*(EXP(1./(AA*tdim(i,j)**2+ &BB*tdim(i,j)+CC))))/etaref RA(I,J) = 1e6 + 0.9e6 pr(i,j) = Pr0 502 CONTINUE DO 530 I=2,L2 DO 530 J=2,M2 C COMPONENTE U ********************** IF(NF.EQ.1) THEN if(partes(i,j).eq.1)gam(i,j) = 1E30 if(partes(i,j).eq.0)gam(i,j) = etai(i,j) CON(I,J) = UANT(I,J)/DTPO AP(I,J) = -1./DTPO ENDIF C COMPONENTE V ********************** IF(NF.EQ.2) THEN if(partes(i,j).eq.1)gam(i,j) = 1E30 if(partes(i,j).eq.0)gam(i,j) = etai(i,j) IF (i.EQ.1) TM = T(I,J) IF (i.NE.1) TM = 0.5*(T(I,J)+T(I,J-1)) IF (j.GT.2) CON(I,J) = (VANT(I,J)/DTpo) & + sqrt(Ra0/Pr0)*(TM) AP(I,J) = -1./DTPO ENDIF 530 CONTINUE C**** TEMPERATURA ********************** DO 540 I=1,L1 DO 540 J=1,M1 IF(NF.EQ.4)THEN rho(i,j) = SQRT(RA(I,J)*PR(I,J)) if(partes(i,j).eq.1)then !solido gam(i,j) = 1. CON(I,J) = TANT(I,J)*rho(i,j)*(1./0.83)/DTPO AP(I,J) = - rho(i,j)*(1./0.83)/DTPO endif if(partes(i,j).eq.0)then !liquido gam(i,j) = 1. CON(I,J) = TANT(I,J)*rho(i,j)/DTPO AP(I,J) = - rho(i,j)/DTPO endif ENDIF 540 CONTINUE C**** CONCENTRACION ********************** DO 550 I=1,L1 DO 550 J=1,M1 IF(NF.EQ.5)THEN rhocp(i,j) = SQRT(RA(I,J)*PR(I,J)) if(partes(i,j).eq.1)then gam(i,j) = 1. CON(I,J) = cANT(I,J)*rhocp(i,j)*(1./0.83)/DTPO AP(I,J) = - rhocp(i,j)*(1./0.83)/DTPO endif if(partes(i,j).eq.0)then gam(i,j) = 1. CON(I,J) = cANT(I,J)*rhocp(i,j)/DTPO AP(I,J) = - rhocp(i,j)/DTPO endif ENDIF 550 CONTINUE RETURN END C*********************************************************************** C----------------METODO DE VOLÚMENES FINITOS 2D C-ACOPLAMIENTO DE ECUACIONES DE MOMENTUM CON EL ALGORITMO SIMPLE-------- C C_______________________ VOLUMENES DE CONTROL--------------------------- C*********************************************************************** LOGICAL LSTOP LINE0110 COMMON/CNTL/LSTOP LINE0120 CALL OPENF CALL GRID LINE0130 CALL SETUP1 LINE0140 CALL START LINE0150 10 CALL DENSE LINE0160 CALL BOUND LINE0170 CALL OUTPUT LINE0180 IF(LSTOP) THEN STOP END IF LINE0190 CALL SETUP2 LINE0200 GO TO 10 END C*********************************************************************** C*********************************************************************** SUBROUTINE DIFLOW LINE0230 COMMON/COEF/FLOW,DIFF,ACOF LINE0240 ACOF=DIFF LINE0250 IF(FLOW.EQ.0.) RETURN LINE0260 TEMP=DIFF-ABS(FLOW)*0.1 LINE0270 ACOF=0. LINE0280 IF(TEMP.LE.0.)RETURN LINE0290 TEMP=TEMP/DIFF LINE0300 ACOF=DIFF*TEMP**5 LINE0310 RETURN LINE0320 END LINE0330 SUBROUTINE SOLVE LINE0340 CHARACTER*10 TITLE LINE0350 LOGICAL LSOLVE,LPRINT,LBLK LINE0360 PARAMETER (NOD=180) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), LINE0365 1 CON(NOD,NOD), LINE0370 2 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), LINE0380 3 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), LINE0390 4 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), LINE0400 5 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), LINE0410 6 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), LINE0420 7 YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), LINE0430 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) LINE0440 COMMON/CHAR/TITLE(13) LINE0450 COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, LINE0460 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, LINE0470 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON LINE0480 ISTF=IST-1 LINE0490 JSTF=JST-1 LINE0500 IT1=L2+IST LINE0510 IT2=L3+IST LINE0520 JT1=M2+JST LINE0530 JT2=M3+JST LINE0540 DO 999 NT=1,NTIMES(NF) LINE0550 DO 999 N=NF,NF LINE0560 C---------------------------------------------------------------------- LINE0570 IF(.NOT.LBLK(NF)) GO TO 10 LINE0580 PT(ISTF)=0. LINE0590 QT(ISTF)=0. LINE0600 DO 11 I=IST,L2 LINE0610 BL=0. LINE0620 BLP=0. LINE0630 BLM=0. LINE0640 BLC=0. LINE0650 DO 12 J=JST,M2 LINE0660 BL=BL+AP(I,J) LINE0670 IF(J.NE.M2) BL=BL-AJP(I,J) LINE0680 IF(J.NE.JST) BL=BL-AJM(I,J) LINE0690 BLP=BLP+AIP(I,J) LINE0700 BLM=BLM+AIM(I,J) LINE0710 BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE0720 1 +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N) LINE0730 12 CONTINUE LINE0740 DENOM=BL-PT(I-1)*BLM LINE0750 IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E20 LINE0760 PT(I)=BLP/DENOM LINE0770 QT(I)=(BLC+BLM*QT(I-1))/DENOM LINE0780 11 CONTINUE LINE0790 BL=0. LINE0800 DO 13 II=IST,L2 LINE0810 I=IT1-II LINE0820 BL=BL*PT(I)+QT(I) LINE0830 DO 13 J=JST,M2 LINE0840 13 F(I,J,N)=F(I,J,N)+BL LINE0850 PT(JSTF)=0. LINE0860 QT(JSTF)=0. LINE0870 DO 21 J=JST,M2 LINE0880 BL=0. LINE0890 BLP=0. LINE0900 BLM=0. LINE0910 BLC=0. LINE0920 DO 22 I=IST,L2 LINE0930 BL=BL+AP(I,J) LINE0940 IF(I.NE.L2) BL=BL-AIP(I,J) LINE0950 IF(I.NE.IST) BL=BL-AIM(I,J) LINE0960 BLP=BLP+AJP(I,J) LINE0970 BLM=BLM+AJM(I,J) LINE0980 BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE0990 1 +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N) LINE1000 22 CONTINUE LINE1010 DENOM=BL-PT(J-1)*BLM LINE1020 IF(ABS(DENOM/BL).LT.1.E-10) DENOM=1.E20 LINE1030 PT(J)=BLP/DENOM LINE1040 QT(J)=(BLC+BLM*QT(J-1))/DENOM LINE1050 21 CONTINUE LINE1060 BL=0. LINE1070 DO 23 JJ=JST,M2 LINE1080 J=JT1-JJ LINE1090 BL=BL*PT(J)+QT(J) LINE1100 DO 23 I=IST,L2 LINE1110 23 F(I,J,N)=F(I,J,N)+BL LINE1120 10 CONTINUE LINE1130 DO 90 J=JST,M2 LINE1140 PT(ISTF)=0. LINE1150 QT(ISTF)=F(ISTF,J,N) LINE1160 DO 70 I=IST,L2 LINE1170 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) LINE1180 PT(I)=AIP(I,J)/DENOM LINE1190 TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N) LINE1200 QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM LINE1210 70 CONTINUE LINE1220 DO 80 II=IST,L2 LINE1230 I=IT1-II LINE1240 80 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) LINE1250 90 CONTINUE LINE1260 DO 190 JJ=JST,M3 LINE1270 J=JT2-JJ LINE1280 PT(ISTF)=0. LINE1290 QT(ISTF)=F(ISTF,J,N) LINE1300 DO 170 I=IST,L2 LINE1310 DENOM=AP(I,J)-PT(I-1)*AIM(I,J) LINE1320 PT(I)=AIP(I,J)/DENOM LINE1330 TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N) LINE1340 QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM LINE1350 170 CONTINUE LINE1360 DO 180 II=IST,L2 LINE1370 I=IT1-II LINE1380 180 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I) LINE1390 190 CONTINUE LINE1400 DO 290 I=IST,L2 LINE1410 PT(JSTF)=0. LINE1420 QT(JSTF)=F(I,JSTF,N) LINE1430 DO 270 J=JST,M2 LINE1440 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) LINE1450 PT(J)=AJP(I,J)/DENOM LINE1460 TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE1470 QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM LINE1480 270 CONTINUE LINE1490 DO 280 JJ=JST,M2 LINE1500 J=JT1-JJ LINE1510 280 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) LINE1520 290 CONTINUE LINE1530 DO 390 II=IST,L3 LINE1540 I=IT2-II LINE1550 PT(JSTF)=0. LINE1560 QT(JSTF)=F(I,JSTF,N) LINE1570 DO 370 J=JST,M2 LINE1580 DENOM=AP(I,J)-PT(J-1)*AJM(I,J) LINE1590 PT(J)=AJP(I,J)/DENOM LINE1600 TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N) LINE1610 QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM LINE1620 370 CONTINUE LINE1630 DO 380 JJ=JST,M2 LINE1640 J=JT1-JJ LINE1650 380 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J) LINE1660 390 CONTINUE LINE1670 999 CONTINUE LINE1680 DO 400 J=2,M2 LINE1690 DO 400 I=2,L2 LINE1700 CON(I,J)=0. LINE1710 AP(I,J)=0. LINE1720 400 CONTINUE LINE1730 RETURN LINE1740 END LINE1750 SUBROUTINE SETUP LINE1760 CHARACTER*10 TITLE LINE1770 LOGICAL LSOLVE,LPRINT,LBLK,LSTOP LINE1780 PARAMETER (NOD=180) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), LINE1785 1 CON(NOD,NOD), LINE1790 2 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), LINE1800 3 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), LINE1810 4 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), LINE1820 5 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), LINE1830 6 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), LINE1840 7 YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), LINE1850 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) LINE1860 COMMON/CHAR/TITLE(13) LINE1870 COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, LINE1880 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, LINE1890 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON LINE1900 COMMON/CNTL/LSTOP LINE1910 COMMON/SORC/SMAX,SSUM LINE1920 COMMON/COEF/FLOW,DIFF,ACOF LINE1930 REAL U(NOD,NOD),V(NOD,NOD),PC(NOD,NOD) LINE1940 EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) LINE1950 101 FORMAT(//,3X,'COMPUTATION IN CARTESIAN COORDINATES') LINE1960 2 FORMAT(//,3X,'COMPUTATION FOR AXISYMMETRIC SITUATION') LINE1970 3 FORMAT(//,3X,'COMPUTATION IN POLAR COORDINATES') LINE1980 4 FORMAT(2X,40(1H*),//) SI01990 ENTRY SETUP1 LINE2000 L2=L1-1 LINE2010 L3=L2-1 LINE2020 M2=M1-1 LINE2030 M3=M2-1 LINE2040 X(1)=XU(2) LINE2050 DO 5 I=2,L2 LINE2060 5 X(I)=0.5*(XU(I+1)+XU(I)) LINE2070 X(L1)=XU(L1) LINE2080 Y(1)=YV(2) LINE2090 DO 10 J=2,M2 LINE2100 10 Y(J)=0.5*(YV(J+1)+YV(J)) LINE2110 Y(M1)=YV(M1) LINE2120 DO 15 I=2,L1 LINE2130 15 XDIF(I)=X(I)-X(I-1) LINE2140 DO 18 I=2,L2 LINE2150 18 XCV(I)=XU(I+1)-XU(I) LINE2160 DO 20 I=3,L2 LINE2170 20 XCVS(I)=XDIF(I) LINE2180 XCVS(3)=XCVS(3)+XDIF(2) LINE2190 XCVS(L2)=XCVS(L2)+XDIF(L1) LINE2200 DO 22 I=3,L3 LINE2210 XCVI(I)=0.5*XCV(I) LINE2220 22 XCVIP(I)=XCVI(I) LINE2230 XCVIP(2)=XCV(2) LINE2240 XCVI(L2)=XCV(L2) LINE2250 DO 35 J=2,M1 LINE2260 35 YDIF(J)=Y(J)-Y(J-1) LINE2270 DO 40 J=2,M2 LINE2280 40 YCV(J)=YV(J+1)-YV(J) LINE2290 DO 45 J=3,M2 LINE2300 45 YCVS(J)=YDIF(J) LINE2310 YCVS(3)=YCVS(3)+YDIF(2) LINE2320 YCVS(M2)=YCVS(M2)+YDIF(M1) LINE2330 IF(MODE.NE.1) GO TO 55 LINE2340 DO 52 J=1,M1 LINE2350 RMN(J)=1.0 LINE2360 52 R(J)=1.0 LINE2370 GO TO 56 LINE2380 55 DO 50 J=2,M1 LINE2390 50 R(J)=R(J-1)+YDIF(J) LINE2400 RMN(2)=R(1) LINE2410 DO 60 J=3,M2 LINE2420 60 RMN(J)=RMN(J-1)+YCV(J-1) LINE2430 RMN(M1)=R(M1) LINE2440 56 CONTINUE LINE2450 DO 57 J=1,M1 LINE2460 SX(J)=1. LINE2470 SXMN(J)=1. LINE2480 IF(MODE.NE.3) GO TO 57 LINE2490 SX(J)=R(J) LINE2500 IF(J.NE.1) SXMN(J)=RMN(J) LINE2510 57 CONTINUE LINE2520 DO 62 J=2,M2 LINE2530 YCVR(J)=R(J)*YCV(J) LINE2540 ARX(J)=YCVR(J) LINE2550 IF(MODE.NE.3) GO TO 62 LINE2560 ARX(J)=YCV(J) LINE2570 62 CONTINUE LINE2580 DO 64 J=4,M3 LINE2590 64 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J) LINE2600 YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3) LINE2610 YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2) LINE2620 IF(MODE.NE.2) GO TO 67 LINE2630 DO 65 J=3,M3 LINE2640 ARXJ(J)=0.25*(1.+RMN(J)/R(J))*ARX(J) LINE2650 65 ARXJP(J)=ARX(J)-ARXJ(J) LINE2660 GO TO 68 LINE2670 67 DO 66 J=3,M3 LINE2680 ARXJ(J)=0.5*ARX(J) LINE2690 66 ARXJP(J)=ARXJ(J) LINE2700 68 ARXJP(2)=ARX(2) LINE2710 ARXJ(M2)=ARX(M2) LINE2720 DO 70 J=3,M3 LINE2730 FV(J)=ARXJP(J)/ARX(J) LINE2740 70 FVP(J)=1.-FV(J) LINE2750 DO 85 I=3,L2 LINE2760 FX(I)=0.5*XCV(I-1)/XDIF(I) LINE2770 85 FXM(I)=1.-FX(I) LINE2780 FX(2)=0. LINE2790 FXM(2)=1. LINE2800 FX(L1)=1. LINE2810 FXM(L1)=0. LINE2820 DO 90 J=3,M2 LINE2830 FY(J)=0.5*YCV(J-1)/YDIF(J) LINE2840 90 FYM(J)=1.-FY(J) LINE2850 FY(2)=0. LINE2860 FYM(2)=1. LINE2870 FY(M1)=1. LINE2880 FYM(M1)=0. LINE2890 CCON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE LINE2900 DO 95 J=1,M1 LINE2910 DO 95 I=1,L1 LINE2920 PC(I,J)=0. LINE2930 U(I,J)=0. LINE2940 V(I,J)=0. LINE2950 CON(I,J)=0. LINE2960 AP(I,J)=0. LINE2970 RHO(I,J)=RHOCON LINE2980 P(I,J)=0. LINE2990 95 CONTINUE LINE3000 IF (MODE.EQ.1) then write(*,101) write(1,*)' CALCULANDO EN COORDENADAS CARTESIANAS' end if IF (MODE.EQ.2) then write(*,101) write(1,*)' CALCULANDO EN COORDENADAS CILINDRICAS' end if IF (MODE.EQ.3) then write(*,101) write(1,*)' COORDENADAS ESFERICAS' end if RETURN LINE3050 ENTRY SETUP2 LINE3060 CCOEFFICIENTS FOR THE U EQUATION LINE3070 NF=1 LINE3080 IF(.NOT.LSOLVE(NF)) GO TO 100 LINE3090 IST=3 LINE3100 JST=2 LINE3110 CALL GAMSOR LINE3120 REL=1.-RELAX(NF) LINE3130 DO 102 I=3,L2 LINE3140 FL=XCVI(I)*V(I,2)*RHO(I,1) LINE3150 FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1) LINE3160 FLOW=R(1)*(FL+FLM) LINE3170 DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2) LINE3180 CALL DIFLOW LINE3190 102 AJM(I,2)=ACOF+AMAX1(0.,FLOW) LINE3200 DO 103 J=2,M2 LINE3210 FLOW=ARX(J)*U(2,J)*RHO(1,J) LINE3220 DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J)) LINE3230 CALL DIFLOW LINE3240 AIM(3,J)=ACOF+AMAX1(0.,FLOW) LINE3250 DO 103 I=3,L2 LINE3260 IF(I.EQ.L2) GO TO 104 LINE3270 FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J)) LINE3280 FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE3290 FLOW=ARX(J)*0.5*(FL+FLP) LINE3300 DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J)) LINE3310 GO TO 105 LINE3320 104 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) LINE3330 DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J)) LINE3340 105 CALL DIFLOW LINE3350 AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) LINE3360 AIP(I,J)=AIM(I+1,J)-FLOW LINE3370 IF(J.EQ.M2) GO TO 106 LINE3380 FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) LINE3390 FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)* LINE3400 1 RHO(I-1,J)) LINE3410 GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+ LINE3420 1 1.0E-30)*XCVI(I) LINE3430 GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)* LINE3440 1 GAM(I-1,J)+1.E-30)*XCVIP(I-1) LINE3450 DIFF=RMN(J+1)*2.*(GM+GMM) LINE3460 GO TO 107 LINE3470 106 FL=XCVI(I)*V(I,M1)*RHO(I,M1) LINE3480 FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1) LINE3490 DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1) LINE3500 107 FLOW=RMN(J+1)*(FL+FLM) LINE3510 CALL DIFLOW LINE3520 AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) LINE3530 AJP(I,J)=AJM(I,J+1)-FLOW LINE3540 VOL=YCVR(J)*XCVS(I) LINE3550 APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1)) LINE3560 1/(XCVS(I)*DT) LINE3570 AP(I,J)=AP(I,J)-APT LINE3580 CON(I,J)=CON(I,J)+APT*U(I,J) LINE3590 AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) LINE3600 1/RELAX(NF) LINE3610 CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J) LINE3620 DU(I,J)=VOL/(XDIF(I)*SX(J)) LINE3630 CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J)) LINE3640 DU(I,J)=DU(I,J)/AP(I,J) LINE3650 103 CONTINUE LINE3660 CALL SOLVE LINE3670 100 CONTINUE LINE3680 CCOEFFICIENTS FOR THE V EQUATION------------------ LINE3690 NF=2 LINE3700 IF(.NOT.LSOLVE(NF)) GO TO 200 LINE3710 IST=2 LINE3720 JST=3 LINE3730 CALL GAMSOR LINE3740 REL=1.-RELAX(NF) LINE3750 DO 202 I=2,L2 LINE3760 AREA=R(1)*XCV(I) LINE3770 FLOW=AREA*V(I,2)*RHO(I,1) LINE3780 DIFF=AREA*GAM(I,1)/YCV(2) LINE3790 CALL DIFLOW LINE3800 202 AJM(I,3)=ACOF+AMAX1(0.,FLOW) LINE3810 DO 203 J=3,M2 LINE3820 FL=ARXJ(J)*U(2,J)*RHO(1,J) LINE3830 FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1) LINE3840 FLOW=FL+FLM LINE3850 DIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J)) LINE3860 CALL DIFLOW LINE3870 AIM(2,J)=ACOF+AMAX1(0.,FLOW) LINE3880 DO 203 I=2,L2 LINE3890 IF(I.EQ.L2) GO TO 204 LINE3900 FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE3910 FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)* LINE3920 1 RHO(I,J-1)) LINE3930 GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+ LINE3940 1 1.E-30)*ARXJ(J) LINE3950 GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)* LINE3960 1 GAM(I,J-1)+1.0E-30)*ARXJP(J-1) LINE3970 DIFF=2.*(GM+GMM)/SXMN(J) LINE3980 GO TO 205 LINE3990 204 FL=ARXJ(J)*U(L1,J)*RHO(L1,J) LINE4000 FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1) LINE4010 DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J))LINE4020 205 FLOW=FL+FLM LINE4030 CALL DIFLOW LINE4040 AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) LINE4050 AIP(I,J)=AIM(I+1,J)-FLOW LINE4060 IF(J.EQ.M2) GO TO 206 LINE4070 AREA=R(J)*XCV(I) LINE4080 FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J) LINE4090 FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))*RMN(J+1) LINE4100 FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I) LINE4110 DIFF=AREA*GAM(I,J)/YCV(J) LINE4120 GO TO 207 LINE4130 206 AREA=R(M1)*XCV(I) LINE4140 FLOW=AREA*V(I,M1)*RHO(I,M1) LINE4150 DIFF=AREA*GAM(I,M1)/YCV(M2) LINE4160 207 CALL DIFLOW LINE4170 AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) LINE4180 AJP(I,J)=AJM(I,J+1)-FLOW LINE4190 VOL=YCVRS(J)*XCV(I) LINE4200 SXT=SX(J) LINE4210 IF(J.EQ.M2) SXT=SX(M1) LINE4220 SXB=SX(J-1) LINE4230 IF(J.EQ.3) SXB=SX(1) LINE4240 APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)* LINE4250 10.5*(SXB+SXMN(J)))/(YCVRS(J)*DT) LINE4260 AP(I,J)=AP(I,J)-APT LINE4270 CON(I,J)=CON(I,J)+APT*V(I,J) LINE4280 AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) LINE4290 1/RELAX(NF) LINE4300 CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J) LINE4310 DV(I,J)=VOL/YDIF(J) LINE4320 CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J)) LINE4330 DV(I,J)=DV(I,J)/AP(I,J) LINE4340 203 CONTINUE LINE4350 CALL SOLVE LINE4360 200 CONTINUE LINE4370 CCOEFFICIENTS FOR THE PRESSURE CORRECTION EQUATION--------------------- LINE4380 NF=3 LINE4390 IF(.NOT.LSOLVE(NF)) GO TO 500 LINE4400 IST=2 LINE4410 JST=2 LINE4420 CALL GAMSOR LINE4430 SMAX=0. LINE4440 SSUM=0. LINE4450 DO 410 J=2,M2 LINE4460 DO 410 I=2,L2 LINE4470 VOL=YCVR(J)*XCV(I) LINE4480 410 CON(I,J)=CON(I,J)*VOL LINE4490 DO 402 I=2,L2 LINE4500 ARHO=R(1)*XCV(I)*RHO(I,1) LINE4510 CON(I,2)=CON(I,2)+ARHO*V(I,2) LINE4520 402 AJM(I,2)=0. LINE4530 DO 403 J=2,M2 LINE4540 ARHO=ARX(J)*RHO(1,J) LINE4550 CON(2,J)=CON(2,J)+ARHO*U(2,J) LINE4560 AIM(2,J)=0. LINE4570 DO 403 I=2,L2 LINE4580 IF(I.EQ.L2) GO TO 404 LINE4590 ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE4600 FLOW=ARHO*U(I+1,J) LINE4610 CON(I,J)=CON(I,J)-FLOW LINE4620 CON(I+1,J)=CON(I+1,J)+FLOW LINE4630 AIP(I,J)=ARHO*DU(I+1,J) LINE4640 AIM(I+1,J)=AIP(I,J) LINE4650 GO TO 405 LINE4660 404 ARHO=ARX(J)*RHO(L1,J) LINE4670 CON(I,J)=CON(I,J)-ARHO*U(L1,J) LINE4680 AIP(I,J)=0. LINE4690 405 IF(J.EQ.M2) GO TO 406 LINE4700 ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) LINE4710 FLOW=ARHO*V(I,J+1) LINE4720 CON(I,J)=CON(I,J)-FLOW LINE4730 CON(I,J+1)=CON(I,J+1)+FLOW LINE4740 AJP(I,J)=ARHO*DV(I,J+1) LINE4750 AJM(I,J+1)=AJP(I,J) LINE4760 GO TO 407 LINE4770 406 ARHO=RMN(M1)*XCV(I)*RHO(I,M1) LINE4780 CON(I,J)=CON(I,J)-ARHO*V(I,M1) LINE4790 AJP(I,J)=0. LINE4800 407 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J) LINE4810 PC(I,J)=0. LINE4820 SMAX=AMAX1(SMAX,ABS(CON(I,J))) LINE4830 SSUM=SSUM+CON(I,J) LINE4840 403 CONTINUE LINE4850 CALL SOLVE LINE4860 CCOME HERE TO CORRECT THE PRESSURE AND VELOCITIES--------------------- LINE4870 DO 501 J=2,M2 LINE4880 DO 501 I=2,L2 LINE4890 P(I,J)=P(I,J)+PC(I,J)*RELAX(NP) LINE4900 IF(I.NE.2) U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J)) LINE4910 IF(J.NE.2) V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J)) LINE4920 501 CONTINUE LINE4930 500 CONTINUE LINE4940 CCOEFFICIENTS FOR OTHER EQUATIONS-------------------------------- LINE4950 IST=2 LINE4960 JST=2 LINE4970 DO 600 N=4,NFMAX LINE4980 NF=N LINE4990 IF(.NOT.LSOLVE(NF)) GO TO 600 LINE5000 CALL GAMSOR LINE5010 REL=1.-RELAX(NF) LINE5020 DO 602 I=2,L2 LINE5030 AREA=R(1)*XCV(I) LINE5040 FLOW=AREA*V(I,2)*RHO(I,1) LINE5050 DIFF=AREA*GAM(I,1)/YDIF(2) LINE5060 CALL DIFLOW LINE5070 602 AJM(I,2)=ACOF+AMAX1(0.,FLOW) LINE5080 DO 603 J=2,M2 LINE5090 FLOW=ARX(J)*U(2,J)*RHO(1,J) LINE5100 DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J)) LINE5110 CALL DIFLOW LINE5120 AIM(2,J)=ACOF+AMAX1(0.,FLOW) LINE5130 DO 603 I=2,L2 LINE5140 IF(I.EQ.L2) GO TO 604 LINE5150 FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J)) LINE5160 DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+ LINE5170 1 XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J)) LINE5180 GO TO 605 LINE5190 604 FLOW=ARX(J)*U(L1,J)*RHO(L1,J) LINE5200 DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J)) LINE5210 605 CALL DIFLOW LINE5220 AIM(I+1,J)=ACOF+AMAX1(0.,FLOW) LINE5230 AIP(I,J)=AIM(I+1,J)-FLOW LINE5240 AREA=RMN(J+1)*XCV(I) LINE5250 IF(J.EQ.M2) GO TO 606 LINE5260 FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) LINE5270 DIFF=AREA*2.*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+ LINE5280 1 YCV(J+1)*GAM(I,J)+1.0E-30) LINE5290 GO TO 607 LINE5300 606 FLOW=AREA*V(I,M1)*RHO(I,M1) LINE5310 DIFF=AREA*GAM(I,M1)/YDIF(M1) LINE5320 607 CALL DIFLOW LINE5330 AJM(I,J+1)=ACOF+AMAX1(0.,FLOW) LINE5340 AJP(I,J)=AJM(I,J+1)-FLOW LINE5350 VOL=YCVR(J)*XCV(I) LINE5360 APT=RHO(I,J)/DT LINE5370 AP(I,J)=AP(I,J)-APT LINE5380 CON(I,J)=CON(I,J)+APT*F(I,J,NF) LINE5390 AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)) LINE5400 1/RELAX(NF) LINE5410 CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*F(I,J,NF) LINE5420 603 CONTINUE LINE5430 CALL SOLVE LINE5440 600 CONTINUE LINE5450 TIME=TIME+DT LINE5460 ITER=ITER+1 LINE5470 IF(ITER.GE.LAST) LSTOP=.TRUE. LINE5480 RETURN LINE5490 END LINE5500 SUBROUTINE SUPPLY LINE5510 CHARACTER*10 TITLE LINE5520 LOGICAL LSOLVE,LPRINT,LBLK LINE5530 PARAMETER (NOD=180) !NUMERO MAXIMO DE NODOS COMMON F(NOD,NOD,10),P(NOD,NOD),RHO(NOD,NOD),GAM(NOD,NOD), LINE5535 1 CON(NOD,NOD), LINE5540 2 AIP(NOD,NOD),AIM(NOD,NOD),AJP(NOD,NOD),AJM(NOD,NOD),AP(NOD,NOD), LINE5550 3 X(NOD),XU(NOD),XDIF(NOD),XCV(NOD),XCVS(NOD), LINE5560 4 Y(NOD),YV(NOD),YDIF(NOD),YCV(NOD),YCVS(NOD), LINE5570 5 YCVR(NOD),YCVRS(NOD),ARX(NOD),ARXJ(NOD),ARXJP(NOD), LINE5580 6 R(NOD),RMN(NOD),SX(NOD),SXMN(NOD),XCVI(NOD),XCVIP(NOD), LINE5590 7 YLS(10),M1S(10),M1ST(10),XLS(10),L1S(10),L1ST(10),NSX,NSY COMMON DU(NOD,NOD),DV(NOD,NOD),FV(NOD),FVP(NOD), LINE5600 1 FX(NOD),FXM(NOD),FY(NOD),FYM(NOD),PT(NOD),QT(NOD) LINE5610 COMMON/CHAR/TITLE(13) LINE5620 COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3, LINE5630 1IST,JST,ITER,LAST,RELAX(13),TIME,DT,XL,YL, LINE5640 2IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON LINE5650 REAL U(NOD,NOD),V(NOD,NOD),PC(NOD,NOD) LINE5660 EQUIVALENCE(F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1)) LINE5670 10 FORMAT(6(1H*),3X,A10,3X,6(1H*),/,8X,12(1H-)) LINE5680 20 FORMAT(4H I =,I6,12I9) LINE5690 30 FORMAT(1HJ) LINE5700 40 FORMAT(1X,I2,3X,(1x,13E9.2)) LINE5710 50 FORMAT(1H ) LINE5720 51 FORMAT(1X,'I =',2X,13(I4,5X)) LINE5730 52 FORMAT(1X,'X =',1x,13E9.2) LINE5740 53 FORMAT(1X,'TH =',1x,13E9.2) LINE5750 54 FORMAT(1X,'J =',2X,13(I4,5X)) LINE5760 55 FORMAT(1X,'Y =',1x,13E9.2) LINE5770 ENTRY UGRID LINE5780 XU(2)=0. LINE5790 DX=XL/FLOAT(L1-2) LINE5800 DO 1 I=3,L1 LINE5810 1 XU(I)=XU(I-1)+DX LINE5820 YV(2)=0. LINE5830 DY=YL/FLOAT(M1-2) LINE5840 DO 2 J=3,M1 LINE5850 2 YV(J)=YV(J-1)+DY LINE5860 RETURN LINE5870 ENTRY NUGRID L1ST(1)=L1S(1) DO 5 I=2,NSX 5 L1ST(I)=L1ST(I-1)+L1S(I) L1=L1ST(NSX) DO 6 I=1,NSX 6 XL=XLS(I)+XL XU(2)=0. I=1 DX=XLS(1)/FLOAT(L1S(1)-1) DO 3 J=3,L1 XU(J)=XU(J-1)+DX IF(J-1.EQ.L1ST(I)) THEN I=I+1 DX=XLS(I)/FLOAT(L1S(I)) IF(I.EQ.NSX) DX=XLS(I)/FLOAT(L1S(I)-1) 3 ENDIF M1ST(1)=M1S(1) DO 7 I=2,NSY 7 M1ST(I)=M1ST(I-1)+M1S(I) M1=M1ST(NSY) DO 8 I=1,NSY 8 YL=YLS(I)+YL YV(2)=0. I=1 DY=YLS(1)/FLOAT(M1S(1)-1) DO 4 J=3,M1 YV(J)=YV(J-1)+DY IF(J-1.EQ.M1ST(I)) THEN I=I+1 DY=YLS(I)/FLOAT(M1S(I)) IF(I.EQ.NSY) DY=YLS(I)/FLOAT(M1S(I)-1) 4 ENDIF RETURN ENTRY PRINT LINE5880 IF(.NOT.LPRINT(3)) GO TO 80 LINE5890 CALCULATE THE STREAM FUNCTION LINE5900 F(2,2,3)=0. LINE5910 DO 82 I=2,L1 LINE5920 IF(I.NE.2) F(I,2,3)=F(I-1,2,3)-RHO(I-1,1)*V(I-1,2) LINE5930 1*R(1)*XCV(I-1) LINE5940 DO 82 J=3,M1 LINE5950 RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1) LINE5960 82 F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1) LINE5970 80 CONTINUE LINE5980 C LINE5990 IF(.NOT.LPRINT(NP)) GO TO 90 LINE6000 C LINE6010 CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATION LINE6020 DO 91 J=2,M2 LINE6030 P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3) LINE6040 91 P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2) LINE6050 DO 92 I=2,L2 LINE6060 P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3) LINE6070 92 P(I,M1)=(P(I,M2)*YCVS(M2)-P(I,M3)*YDIF(M1))/YDIF(M2) LINE6080 P(1,1)=P(2,1)+P(1,2)-P(2,2) LINE6090 P(L1,1)=P(L2,1)+P(L1,2)-P(L2,2) LINE6100 P(1,M1)=P(2,M1)+P(1,M2)-P(2,M2) LINE6110 P(L1,M1)=P(L2,M1)+P(L1,M2)-P(L2,M2) LINE6120 PREF=P(IPREF,JPREF) LINE6130 DO 93 J=1,M1 LINE6140 DO 93 I=1,L1 LINE6150 93 P(I,J)=P(I,J)-PREF LINE6160 90 CONTINUE LINE6170 IEND=0 LINE6200 301 IF(IEND.EQ.L1) GO TO 310 LINE6210 IBEG=IEND+1 LINE6220 IEND=IEND+13 LINE6230 IEND=MIN0(IEND,L1) LINE6240 IF(MODE.EQ.3)GO TO 302 LINE6270 GO TO 303 LINE6290 302 continue 303 GO TO 301 LINE6310 310 JEND=0 LINE6320 311 IF(JEND.EQ.M1) GO TO 320 LINE6340 JBEG=JEND+1 LINE6350 JEND=JEND+13 LINE6360 JEND=MIN0(JEND,M1) LINE6370 GO TO 311 LINE6410 320 CONTINUE LINE6420 C LINE6430 DO 999 N=1,NGAM LINE6440 NF=N LINE6450 IF(.NOT.LPRINT(NF))GO TO 999 LINE6460 IFST=1 LINE6490 JFST=1 LINE6500 IF(NF.EQ.1.OR.NF.EQ.3) IFST=2 LINE6510 IF(NF.EQ.2.OR.NF.EQ.3) JFST=2 LINE6520 IBEG=IFST-13 LINE6530 110 CONTINUE LINE6540 IBEG=IBEG+13 LINE6550 IEND=IBEG+12 LINE6560 IEND=MIN0(IEND,L1) LINE6570 JFL=JFST+M1 LINE6610 DO 115 JJ=JFST,M1 LINE6620 J=JFL-JJ LINE6630 115 CONTINUE LINE6650 IF(IEND.LT.L1) GO TO 110 LINE6660 999 CONTINUE LINE6670 RETURN LINE6680 END LINE6690