ckon3
C CKON3 SOURCE CB215821 20/11/25 13:19:49 10792 & IROF,IVITF,IPF,IFRMAF,ISCAF,PROPHY, & ICHPSU,ICHPDI, & MELEMC,MELEMF,MELEFE, & IZG1,IZG2,IZG3,IZG4,IZG5,DT,DIAMEL,NLCEMI, & LOGNC,LOGAN,MESERR) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : CKON3 C C DESCRIPTION : Voir CKON C C Cas deux dimensions, gaz "thermally perfect" C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF C C************************************************************************ C C C APPELES (Outils C CASTEM) : KRIPAD, LICHT C C APPELES (Calcul) : FVLHTE, FCOLTE C C C************************************************************************ C C ENTREES C C C 1) PARAMETRES C C LOGME : (LOGICAL); .TRUE. -> MULTI-ESPECES C .FALSE. -> MONO-ESPECE C C INDMET : 1 VLH (van Leer Hanel FVS) C C 2 Colella-Glaz FDS C C NORDP1 : (ordre des polynoms cv(T)) + 1 C C 2) Pointeurs des MCHAMLs C C IROF : MCHAML sur "FACEL" contenant la masse volumique C ("gauche" et "droite"); C C IVITF : MCHAML sur "FACEL" contenant la vitesse dans le repaire C local (n,t) et les cosinus directeurs des repaire local; C C IPF : MCHAML sur "FACEL" contenant la pression; C C IFRAMAF : MCHAML sur "FACEL", contenant les fractions massiques C si LOGME = .TRUE.; C LOGME = .FALSE. -> IFRAMAF = 0 C C ISCAF : MCHAML sur "FACEL", contenant les scalires passifs a C transporter (ou ISCAF = 0) C C 3) Pointeur sur la table des proprietes de gaz C C PROPHY C C 4) Pointeurs de CHPOINTs de la table DOMAINE C C ICHPSU : CHPOINT "FACE" contenant la surface des faces C C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum C de chaque element C C C 5) Pointeurs de MELEME de la table DOMAINE C C MELEMC : MELEME 'CENTRE' du SPG des CENTRES C C MELEMF : MELEME 'FACE' du SPG des FACES C C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM C C SORTIES (il faudrait dire E/S) C C IZGi : pointeurs de CHPOINTs "FACE" dont les valeurs C se trouvent dans la table KIZX . 'EQEX' . 'KIZG' C aux indices 'IZGi' C C IZG1 : Increment de masse voluique C C IZG2 : Increment de quantite de mouvement C C IZG3 : Increment de l'energie totale C C IZG4 : Increment de les Masse Volumiques des Especies C (si LOGME = .TRUE.) C C IZG5 : Increment de scalaires passifs C C DIAMEL : 'minimum' diametre du maillage C C NLCEMI : numero local du CENTRE ou le diametre est 'minimum' C C DT : pas de temps pour le respect de la CFL-like condition C DT < DIAMEL /2 /max(Lambda_i) C En maillage regulier cette condition garantie la C non-interaction des ondes C C C LOGNC : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée C dans pour la solution du probleme Riemann n'a pas bien C marchéee C C LOGAN : (LOGICAL): si .TRUE. une anomalie à été detectée C C MESERR : pour l'ecriture des messages d'erreurs C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : 22.12.98 Creation C C 08.02.00 Ajout du flux numerique de Colella-Glaz C C 21.02.00 Ajout de transport de scalaires passifs C C************************************************************************ C C C N.B.: On suppose qu'on a déjà controllé RO, P > 0 C GAMMA \in (1,3) C Y \in (0,1) C Si non il faut le faire!!! C C************************************************************************ C C C**** Variables de COOPTIO C C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES C & ,IECHO, IIMPI, IOSPI C & ,IDIM CC & ,MCOORD C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU C & ,NORINC,NORVAL,NORIND,NORVAD C & ,NUCROU, IPSAUV CC IMPLICIT INTEGER(I-N) INTEGER I1, I2 & ,INDMET, NORDP1 & ,IROF,IVITF,IPF,ISCAF,IFRMAF & ,ICHPSU,ICHPDI,MELEMC,MELEMF,MELEFE & ,IGEOMC,IGEOMF & ,IZG1,IZG2,IZG3,IZG4,IZG5,NLCEMI & ,NESP, NFAC, NSCA, NMA, NMA2 & ,NLCF, NGCEG, NGCED, NLCEG, NLCED & ,NGCF, NLCF1, SPG1, SPG2 REAL*8 DIAMEL, DT, UNSDT, CELLT & , ROG, UNG, UTG, PG, GAMG, TG & , ROD, UND, UTD, PD, GAMD, TD & , SURF,CNX, CNY, CTX , CTY & , CELL, DIAMG, DIAMD, DIAM & , ASON, LAMBDA & , Y1G, Y1D, RTOTG, RTOTD, CVTOTG, CVTOTD, ETHERG, ETHERD & , YG, YD, FUNTG, FUNTD, DCV, XST LOGICAL LOGME, LOGNC, LOGAN CHARACTER*(40) MESERR CHARACTER*(8) TYPE C C**** Segment des proprietes du gaz C SEGMENT PROPHY REAL*8 ACV(NORDP1,NESP+1), R(NESP+1), H0K(NESP+1) & ,ACVTOG(NORDP1), ACVTOD(NORDP1) ENDSEGMENT C C**** Les Includes C C C**** LES INCLUDES C -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMCHAML POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, & MELT1X.MELVAL, MELT1Y.MELVAL POINTEUR MELVUN.MELVAL, MELVUT.MELVAL POINTEUR MELRO.MELVAL, MELP.MELVAL -INC SMCHPOI POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL & , MPOVG1.MPOVAL, MPOVG2.MPOVAL & , MPOVG3.MPOVAL, MPOVG4.MPOVAL, MPOVG5.MPOVAL POINTEUR MCHAMY.MCHAML, MCHAMS.MCHAML -INC SMELEME -INC SMLMOTS -INC SMLENTI C C**** Les fractiones massiques. C SEGMENT FRAMAS REAL*8 YET(NMA) ENDSEGMENT POINTEUR FRAMAG.FRAMAS, FRAMAD.FRAMAS, SCAG.FRAMAS, SCAD.FRAMAS C C**** Les flux aux interface dans le repaire (n,t) C SEGMENT IFLUX ENDSEGMENT POINTEUR IFLU1.IFLUX, IFLU2.IFLUX C C C**** Initialisation des MCHAMLs C C**** Masse volumique C MCHEL1 = IROF SEGACT MCHEL1 MCHAM1 = MCHEL1.ICHAML(1) SEGACT MCHAM1 MELRO = MCHAM1.IELVAL(1) SEGDES MCHEL1 SEGDES MCHAM1 C C**** Pression C MCHEL1 = IPF SEGACT MCHEL1 MCHAM1 = MCHEL1.ICHAML(1) SEGACT MCHAM1 MELP = MCHAM1.IELVAL(1) SEGDES MCHEL1 SEGDES MCHAM1 C C**** Vitesse et cosinus directeurs du repere (n,t) C MCHEL1 = IVITF SEGACT MCHEL1 C C**** La vitesse a comme SPG MELEFE C Le cosinus directeurs ont comme SPG MELEMF C C MCHAM1 -> Cosinus directeurs C MCHAM2 -> Vitesse C SPG1 = MCHEL1.IMACHE(1) SPG2 = MCHEL1.IMACHE(2) IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN MCHAM1 = MCHEL1.ICHAML(1) MCHAM2 = MCHEL1.ICHAML(2) ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN MCHAM1 = MCHEL1.ICHAML(2) MCHAM2 = MCHEL1.ICHAML(1) ELSE LOGAN = .TRUE. GOTO 9999 ENDIF SEGACT MCHAM1 MELVNX = MCHAM1.IELVAL(1) MELVNY = MCHAM1.IELVAL(2) MELT1X = MCHAM1.IELVAL(3) MELT1Y = MCHAM1.IELVAL(4) SEGDES MCHAM1 SEGACT MCHAM2 MELVUN = MCHAM2.IELVAL(1) MELVUT = MCHAM2.IELVAL(2) SEGDES MCHAM2 SEGDES MCHEL1 C C**** Fractions massiques C IF(LOGME)THEN MCHEL1 = IFRMAF SEGACT MCHEL1 MCHAMY = MCHEL1.ICHAML(1) SEGACT MCHAMY C C******* Numero d'especes dans les equations d'Euler C NESP = MCHAMY.IELVAL(/1) DO I1 = 1, NESP MELVA1 = MCHAMY.IELVAL(I1) SEGACT MELVA1 ENDDO NMA = NESP SEGINI FRAMAG SEGINI FRAMAD SEGDES MCHEL1 ELSE C C******* Definition minimale de YET, necessaire pour transmetre YET aux C subroutines FORTRAN qui calculent les flux C NMA = 1 SEGINI FRAMAG SEGINI FRAMAD NESP = 0 ENDIF C C**** Scalaires passifs C IF(ISCAF .GT. 0)THEN MCHEL1 = ISCAF SEGACT MCHEL1 MCHAMS = MCHEL1.ICHAML(1) SEGACT MCHAMS NSCA = MCHAMS.IELVAL(/1) DO I1 = 1, NSCA, 1 MELVA1 = MCHAMS.IELVAL(I1) SEGACT MELVA1 ENDDO NMA = NSCA SEGINI SCAG SEGINI SCAD SEGDES MCHEL1 ELSE C C******* Definition minimale de YET, necessaire pour transmetre YET aux C subroutines FORTRAN qui calculent les flux C NMA = 1 SEGINI SCAG SEGINI SCAD NSCA= 0 ENDIF C C**** Initialisation des MELEMEs C C 'CENTRE', 'FACEL' C IPT2 = MELEFE SEGACT IPT2 NFAC = IPT2.NUM(/2) C C**** KRIPAD pour la correspondance global/local de centre C C C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements C C Si i est le numero global d'un noeud de ICEN, C MLENT1.LECT(i) contient sa position, i.e. C C I = numero global du noeud centre C MLENT1.LECT(i) = numero local du noeud centre C C MLENT1 déjà activé, i.e. C C SEGACT MLENT1 C C C**** KRIPAD pour la correspondance global/local de 'FACE' C C C**** Initialisation de flux C NMA2 = 4 + NESP + NSCA SEGINI IFLU1 SEGINI IFLU2 C C**** IFLU2 = segment de travail en FLUVLH; c'est plus rapide le definir ici C C C**** CHPOINTs de la table DOMAINE C C C**** LICHT active les MPOVALs en *MOD C C i.e. C C SEGACT MPOVSU*MOD C SEGACT MPOVDI*MOD C C C**** Les FLUX aux face C C La densité C C C SEGACT MPOVG1*MOD C C**** Les debits C C C SEGACT MPOVG2*MOD C C**** L'energie totale volumique C C C SEGACT MPOVG3*MOD C C**** Les Fractions Massiques C IF(LOGME)THEN C C SEGACT MPOVG4*MOD C ENDIF C C**** Les scalaires passifs C IF(ISCAF .GT. 0)THEN C C SEGACT MPOVG5*MOD C ENDIF C C**** Activation des MCHAMLs C SEGACT MELRO SEGACT MELP SEGACT MELVUN SEGACT MELVUT SEGACT MELVNX SEGACT MELVNY SEGACT MELT1X SEGACT MELT1Y C C**** Initialisation de 1/DT C UNSDT = 0.0D0 C C**** BOUCLE SUR FACEL pour le calcul du FLUX C DO NLCF = 1, NFAC C C******* NLCF = numero local du centre de facel C NGCF = numero global du centre de facel C NLCF1 = numero local du centre de face C NGCEG = numero global du centre ELT "gauche" C NLCEG = numero local du centre ELT "gauche" C NGCED = numero global du centre ELT "droite" C NLCED = numero local du centre ELT "droite" C NGCEG = IPT2.NUM(1,NLCF) NGCED = IPT2.NUM(3,NLCF) NGCF = IPT2.NUM(2,NLCF) NLCF1 = MLENT2.LECT(NGCF) NLCEG = MLENT1.LECT(NGCEG) NLCED = MLENT1.LECT(NGCED) C C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris. C IF(NLCF .NE. NLCF1)THEN MESERR = 'Il ne faut pas jouer avec la console. ' LOGAN = .TRUE. GOTO 9999 ENDIF C C******* Recuperation des Etats "gauche" et "droite" C ROG = MELRO.VELCHE(1,NLCF) UNG = MELVUN.VELCHE(1,NLCF) UTG = MELVUT.VELCHE(1,NLCF) PG = MELP.VELCHE(1,NLCF) C ROD = MELRO.VELCHE(3,NLCF) UND = MELVUN.VELCHE(3,NLCF) UTD = MELVUT.VELCHE(3,NLCF) PD = MELP.VELCHE(3,NLCF) C CNX = MELVNX.VELCHE(1,NLCF) CNY = MELVNY.VELCHE(1,NLCF) CTX = MELT1X.VELCHE(1,NLCF) CTY = MELT1Y.VELCHE(1,NLCF) C C******* Le fractiones massiques, R et ACVTO C RTOTG = 0.0D0 RTOTD = 0.0D0 CVTOTG = 0.0D0 CVTOTD = 0.0D0 Y1G = 1.0D0 Y1D = 1.0D0 ETHERG = 0.0D0 ETHERD = 0.0D0 DO I1= 1, NORDP1, 1 PROPHY.ACVTOG(I1) =0.0D0 PROPHY.ACVTOD(I1) =0.0D0 ENDDO DO I1 = 1, NESP, 1 MELVA1 = MCHAMY.IELVAL(I1) YG = MELVA1.VELCHE(1,NLCF) YD = MELVA1.VELCHE(3,NLCF) Y1G = Y1G - YG Y1D = Y1D - YD FRAMAG.YET(I1) = YG FRAMAD.YET(I1) = YD RTOTG = RTOTG + YG * PROPHY.R(I1) RTOTD = RTOTD + YD * PROPHY.R(I1) ENDDO ENDDO RTOTG = RTOTG + Y1G * PROPHY.R(NESP+1) RTOTD = RTOTD + Y1D * PROPHY.R(NESP+1) ENDDO TG = PG / (ROG * RTOTG) TD = PD / (ROD * RTOTD) FUNTG = 1.0D0 FUNTD = 1.0D0 CVTOTG = PROPHY.ACVTOG(1) ETHERG = CVTOTG * TG CVTOTD = PROPHY.ACVTOD(1) ETHERD = CVTOTD * TD DO I1 = 2, NORDP1, 1 FUNTG = FUNTG * TG FUNTD = FUNTD * TD DCV = PROPHY.ACVTOG(I1) * FUNTG CVTOTG = CVTOTG + DCV ETHERG = ETHERG + DCV * TG / I1 DCV = PROPHY.ACVTOD(I1) * FUNTD CVTOTD = CVTOTD + DCV ETHERD = ETHERD + DCV * TD / I1 ENDDO GAMG = (CVTOTG + RTOTG) /CVTOTG GAMD = (CVTOTD + RTOTD) /CVTOTD C C******* Les scalaires passifs C DO I1 = 1, NSCA, 1 MELVA1 = MCHAMS.IELVAL(I1) SCAG.YET(I1) = MELVA1.VELCHE(1,NLCF) SCAD.YET(I1) = MELVA1.VELCHE(3,NLCF) ENDDO C C******* On a defini (ROg,ROUNg,ROUTg,Pg,(Yg)), (ROd,ROUNd,ROUTd,Pd,(Yd)) C et on a déjà verifié ROg, ROd, Pg, Pd > 0 et 0<Y_i<1 C C C******* Calcul du flux aux interfaces C IF(INDMET .EQ. 1)THEN C C********** VLH FVS C & GAMG,ROG,PG,UNG,UTG,ETHERG, & GAMD,ROD,PD,UND,UTD,ETHERD, & FRAMAG.YET,FRAMAD.YET, & SCAG.YET,SCAD.YET, & CELLT) ELSEIF(INDMET .EQ. 2)THEN C C******* Colella-Glaz FDS (avec Entropy-fix) C XST=0.0D0 & GAMG,RTOTG,PROPHY.ACVTOG,ROG,PG,TG,UNG,UTG,ETHERG, & GAMD,RTOTD,PROPHY.ACVTOD,ROD,PD,TD,UND,UTD,ETHERD, & FRAMAG.YET,FRAMAD.YET, & SCAG.YET,SCAD.YET, & IFLU1.FLUX,CELLT, & LOGNC,MESERR,LOGAN) ENDIF C IF(LOGAN) GOTO 9999 IF(LOGNC) GOTO 9999 C C******* Ecriture des flux C C FLUX(1) = RO Un RO Un C FLUX(2) = RO Un Un + P -> RO Un Ux + P CNX C FLUX(3) = RO Un Ut -> RO Un Uy + P CNY C FLUX(4) = RO Un Et RO Un Et C SURF = MPOVSU.VPOCHA(NLCF,1) MPOVG1.VPOCHA(NLCF,1) = MPOVG1.VPOCHA(NLCF,1) + MPOVG2.VPOCHA(NLCF,1) = MPOVG2.VPOCHA(NLCF,1) + MPOVG2.VPOCHA(NLCF,2) = MPOVG2.VPOCHA(NLCF,2) + MPOVG3.VPOCHA(NLCF,1) = MPOVG3.VPOCHA(NLCF,1) + DO I1 = 1, NESP, 1 & * SURF ENDDO DO I1 = 1, NSCA, 1 & * SURF ENDDO C C******* Calcul du pas du temps (CFL) C C****** a) etat a l'interface C DIAMG = MPOVDI.VPOCHA(NLCEG,1) DIAMD = MPOVDI.VPOCHA(NLCED,1) DIAM = MIN(DIAMG,DIAMD) CELL = 1.0D0/DIAM/CELLT IF(CELL .GT. UNSDT)THEN UNSDT = CELL DIAMEL = DIAMG NLCEMI = NLCEG ENDIF C C****** b) etat gauche C ASON = SQRT(GAMG*PG/ROG) LAMBDA = ABS(UNG) + ASON CELL = LAMBDA / DIAM IF(CELL .GT. UNSDT)THEN UNSDT = CELL DIAMEL = DIAMG NLCEMI = NLCEG ENDIF C C****** C) etat droite C ASON = SQRT(GAMD*PD/ROD) LAMBDA = ABS(UND) + ASON CELL = LAMBDA / DIAM IF(CELL .GT. UNSDT)THEN UNSDT = CELL DIAMEL = DIAMD NLCEMI = NLCED ENDIF C C C**** Fin boucle sur FACEL C ENDDO C C**** Pas du temps (condition de non interaction en 1D) C DT = 0.5D0 / UNSDT C C**** Desactivation des segments et C on detruit les MCHAMLs C C C**** SEGSUP FRAMAG C SEGSUP FRAMAD C C meme si LOGME = .FALSE. C SEGSUP FRAMAG SEGSUP FRAMAD C SEGSUP MLENT1 SEGDES MLENT2 SEGDES IPT2 C SEGSUP IFLU1 SEGSUP IFLU2 C SEGDES MPOVSU SEGDES MPOVDI C SEGDES MPOVG1 SEGDES MPOVG2 SEGDES MPOVG3 C SEGDES MELRO SEGDES MELP SEGDES MELVUN SEGDES MELVUT SEGDES MELVNX SEGDES MELVNY SEGDES MELT1X SEGDES MELT1Y C IF(LOGME) THEN DO I1 = 1, NESP MELVA1 = MCHAMY.IELVAL(I1) SEGDES MELVA1 ENDDO SEGDES MPOVG4 C SEGDES MCHAMY ENDIF IF(ISCAF .GT. 0) THEN DO I1 = 1, NSCA, 1 MELVA1 = MCHAMS.IELVAL(I1) SEGDES MELVA1 ENDDO SEGDES MPOVG5 C SEGDES MCHAMS ENDIF C 9999 CONTINUE C RETURN END C
© Cast3M 2003 - Tous droits réservés.
Mentions légales