C KFMM1 SOURCE OF166741 24/12/13 21:16:09 12097 SUBROUTINE KFMM1(IRN,IGN,IRETN,IGAMN, & ICHPSU,ICHPDI,ICHPVO,INORM, & MELEMC,MELEMF,MELEFE,DCFL, & MELLIM,ICHLIM, & ICHRES) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : KFMM1 C C DESCRIPTION : Voir KFMM C C Cas deux dimensions, gaz "calorically perfect" C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF C C************************************************************************ C ENTREES C C C 1) Pointeurs des CHPOINTs C C IRN : CHPOINT 'CENTRE' contenant la masse volumique C C IGN : CHPOINT 'CENTRE' contenant la q.d.m. C C IRETN : CHPOINT 'CENTRE' contenant l'energie totale C C IGAMN : CHPOINT 'CENTRE' contenant le gamma C C ICHLIM : CHPOINT conditions aux bords C C 3) 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 ICHPVO : CHPOINT "CENTRE" contenant le volume C de chaque element C C INORM : CHPOINT "CENTRE" contenant le normales aux faces C C C 4) 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 MELLIM : MELEME SPG conditions aux bords C C 5) C C DCFL = le double de la CFL C C C SORTIES C C ICHRES : resultat C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : Avril 2002 : creation C Janvier 2003: implementation de condition aux limites 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 IMPLICIT INTEGER(I-N) INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPDI,ICHPVO,INORM & ,MELEMC,MELEMF,MELEFE,ICHRES,ISPG1,ISPG2,NFAC & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED & ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, REC, PG, VOLG, DIAMG & , ROD, RUXD, RUYD, UND, RETD, GAMD, PD, VOLD, DIAMD & , CNX, CNY, DCFL & , SURF, SIGMAF & , UNSDT, UNSDTF, UXD, UYD CHARACTER*8 TYPE C C**** LES INCLUDES C -INC PPARAM -INC CCOPTIO -INC SMCHPOI POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL & , MPORES.MPOVAL, MPODTI.MPOVAL, MPOVOL.MPOVAL & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL & , MPNORM.MPOVAL, MPLIM.MPOVAL -INC SMELEME -INC SMLMOTS -INC SMLENTI POINTEUR MLELIM.MLENTI C C**** Initialisation des MLENTI des conditions aux limites C C CALL KRIPAD(MELLIM,MLELIM) C SEGINI MLELIM 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 CALL KRIPAD(MELEMC,MLENT1) 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 SEGINI MLENT1 C C C**** KRIPAD pour la correspondance global/local de 'FACE' C CALL KRIPAD(MELEMF,MLENT2) C SEGINI MLENT2 C C C**** CHPOINTs de la table DOMAINE C CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF) CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC) CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC) CALL LICHT(INORM,MPNORM,TYPE,IGEOMC) C C**** LICHT active les MPOVALs en *MOD C C i.e. C C SEGACT MPOVSU*MOD C SEGACT MPOVOL*MOD C SEGACT MPOVDI*MOD C SEGACT MPNORM*MOD C CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC) C SEGACT MPORES*MOD C C MPODTI initialisé a zero; MPODTI = 1 / DT C SEGINI, MPODTI=MPORES C CALL LICHT(IRN,MPRN,TYPE,IGEOMC) CALL LICHT(IGN,MPGN,TYPE,IGEOMC) CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC) CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC) CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM) C C SEGACT MPRN*MOD C SEGACT MPGN*MOD C SEGACT MPRETN*MOD C SEGACT MPGAMN*MOD C SEGACT MPLIM*MOD C IF(IGEOMF .NE. MELEMF)THEN WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso' WRITE(IOIMP,*) 'Probleme de SPG' C 21 2 C Données incompatibles CALL ERREUR(21) GOTO 9999 ENDIF IF(IGEOMC .NE. MELEMC)THEN WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso' WRITE(IOIMP,*) 'Probleme de SPG' C 21 2 C Données incompatibles CALL ERREUR(21) GOTO 9999 ENDIF C 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 WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso' WRITE(IOIMP,*) 'Probleme de SPG' C 21 2 C Données incompatibles CALL ERREUR(21) GOTO 9999 ENDIF C CNX = MPNORM.VPOCHA(NLCF,1) CNY = MPNORM.VPOCHA(NLCF,2) C C******* Recuperation des Etats "gauche" et "droite" C C ROG = MPRN.VPOCHA(NLCEG,1) RUXG = MPGN.VPOCHA(NLCEG,1) RUYG = MPGN.VPOCHA(NLCEG,2) UNG = (RUXG * CNX) + (RUYG * CNY) UNG = UNG / ROG RETG = MPRETN.VPOCHA(NLCEG,1) GAMG = MPGAMN.VPOCHA(NLCEG,1) REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG)) REC = REC / ROG PG = (GAMG - 1.0D0) * (RETG - REC) VOLG = MPOVOL.VPOCHA(NLCEG,1) DIAMG = MPOVDI.VPOCHA(NLCEG,1) C ROD = MPRN.VPOCHA(NLCED,1) RUXD = MPGN.VPOCHA(NLCED,1) RUYD = MPGN.VPOCHA(NLCED,2) RETD = MPRETN.VPOCHA(NLCED,1) GAMD = MPGAMN.VPOCHA(NLCED,1) REC = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD)) REC = REC / ROD PD = (GAMD - 1.0D0) * (RETD - REC) VOLD = MPOVOL.VPOCHA(NLCED,1) DIAMD = MPOVDI.VPOCHA(NLCED,1) IF(NLCEG .NE. NLCED)THEN UND = (RUXD * CNX) + (RUYD * CNY) UND = UND / ROD ELSE C Murs au condition aux limite NLLIM=MLELIM.LECT(NGCF) IF(NLLIM .EQ.0)THEN C Mur UND = -1.0D0 * UNG ELSE ROD = MPLIM.VPOCHA(NLLIM,1) UXD = MPLIM.VPOCHA(NLLIM,2) UYD = MPLIM.VPOCHA(NLLIM,3) PD = MPLIM.VPOCHA(NLLIM,4) GAMD = GAMG UND = (UXD * CNX) + (UYD * CNY) RETD = ((1.0D0/(GAMD - 1.0D0))*PD)+ & (0.5D0*ROD*((UXD*UXD)+(UYD*UYD))) VOLD = VOLG ENDIF ENDIF C SURF = MPOVSU.VPOCHA(NLCF,1) SIGMAF = (0.5D0 * (GAMG + GAMD)) * (PG + PD) / (ROG + ROD) SIGMAF = SIGMAF ** 0.5D0 SIGMAF = (0.5D0 * (ABS(UNG) + ABS(UND))) + SIGMAF 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