1 : C KFMM2 SOURCE CHAT 05/01/13 00:55:45 5004 2 : SUBROUTINE KFMM2(IRN,IGN,IRETN,IGAMN, 3 : & ICHPSU,ICHPVO,INORM, 4 : & MELEMC,MELEMF,MELEFE, 5 : & MELLIM,ICHLIM, 6 : & ICHRES) 7 : C************************************************************************ 8 : C 9 : C PROJET : CASTEM 2000 10 : C 11 : C NOM : KFMM2 12 : C 13 : C DESCRIPTION : Voir KFMM 14 : C 14 : C Cas deux dimensions, gaz "calorically perfect" 15 : C 16 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 17 : C 18 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF 19 : C 20 : C************************************************************************ 21 : C ENTREES 22 : C 23 : C 24 : C 25 : C 1) Pointeurs des CHPOINTs 26 : C 27 : C IRN : CHPOINT 'CENTRE' contenant la masse volumique 28 : C 29 : C IGN : CHPOINT 'CENTRE' contenant la q.d.m. 30 : C 31 : C IRETN : CHPOINT 'CENTRE' contenant l'energie totale 32 : C 33 : C IGAMN : CHPOINT 'CENTRE' contenant le gamma 34 : C 35 : C ICHLIM : CHPOINT conditions aux bords 36 : C 37 : C 3) Pointeurs de CHPOINTs de la table DOMAINE 38 : C 39 : C ICHPSU : CHPOINT "FACE" contenant la surface des faces 40 : C 41 : C ICHPVO : CHPOINT "CENTRE" contenant le volume 42 : C de chaque element 43 : C 44 : C 45 : C 4) Pointeurs de MELEME de la table DOMAINE 46 : C 47 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES 48 : C 49 : C MELEMF : MELEME 'FACE' du SPG des FACES 50 : C 51 : C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM 52 : C 53 : C MELLIM : MELEME SPG conditions aux bords 54 : C 55 : C SORTIES 56 : C 57 : C ICHRES : resultat 58 : C 59 : C************************************************************************ 60 : C 61 : C HISTORIQUE (Anomalies et modifications éventuelles) 62 : C 63 : C HISTORIQUE : Avril 2002 : creation 64 : C Janvier 2003: implementation de condition aux limites 65 : C 66 : C************************************************************************ 67 : C 68 : C 69 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0 70 : C GAMMA \in (1,3) 71 : C Y \in (0,1) 72 : C Si non il faut le faire!!! 73 : C 74 : C************************************************************************ 75 : C 76 : C 77 : C**** Variables de COOPTIO 78 : C 79 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 80 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 81 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 82 : C & ,IECHO, IIMPI, IOSPI 83 : C & ,IDIM 84 : C & ,MCOORD 85 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 86 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 87 : C & ,NORINC,NORVAL,NORIND,NORVAD 88 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI 89 : C 90 : IMPLICIT INTEGER(I-N) 91 : INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPVO,INORM 92 : & ,MELEMC,MELEMF,MELEFE,ICHRES,ISPG1,ISPG2,NFAC 93 : & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED 94 : & ,MELLIM,ICHLIM,NLLIM 95 : REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, REC, PG, VOLG, UTG 96 : & , ROD, RUXD, RUYD, UND, RETD, GAMD, PD, VOLD, UTD 97 : & , CNX, CNY, CTX, CTY 98 : & , SURF, SIGMAF 99 : & , UXG, UYG, RHTG 100 : & , UXD, UYD, RHTD, COEF 101 : & ,FR, FRUX, FRUY, FRET 102 : 103 : CHARACTER*8 TYPE 104 : C 105 : C**** LES INCLUDES 106 : C 107 : -INC CCOPTIO 108 : -INC SMCHPOI 109 : POINTEUR MPOVSU.MPOVAL 110 : & , MPORES.MPOVAL, MPOVOL.MPOVAL 111 : & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL 112 : & , MPNORM.MPOVAL, MPLIM.MPOVAL 113 : -INC SMELEME 114 : -INC SMLMOTS 115 : -INC SMLENTI 116 : POINTEUR MLELIM.MLENTI 117 : C 118 : C**** Initialisation des MLENTI des conditions aux limites 119 : C 120 : C 120 : CALL KRIPAD(MELLIM,MLELIM) 121 : C SEGINI MLELIM 122 : C 123 : C**** Initialisation des MELEMEs 124 : C 125 : C 'CENTRE', 'FACEL' 126 : C 127 : IPT2 = MELEFE 128 : SEGACT IPT2 129 : NFAC = IPT2.NUM(/2) 130 : C 131 : C**** KRIPAD pour la correspondance global/local de centre 132 : C 132 : CALL KRIPAD(MELEMC,MLENT1) 133 : C 134 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements 135 : C 136 : C Si i est le numero global d'un noeud de ICEN, 137 : C MLENT1.LECT(i) contient sa position, i.e. 138 : C 139 : C I = numero global du noeud centre 140 : C MLENT1.LECT(i) = numero local du noeud centre 141 : C 142 : C MLENT1 déjà activé, i.e. 143 : C 144 : C SEGINI MLENT1 145 : C 146 : C 147 : C**** KRIPAD pour la correspondance global/local de 'FACE' 148 : C 148 : CALL KRIPAD(MELEMF,MLENT2) 149 : C SEGINI MLENT2 150 : C 151 : C**** CHPOINTs de la table DOMAINE 152 : C 152 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF) 152 : CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC) 152 : CALL LICHT(INORM,MPNORM,TYPE,IGEOMC) 153 : C 154 : C**** LICHT active les MPOVALs en *MOD 155 : C 156 : C i.e. 157 : C 158 : C SEGACT MPOVSU*MOD 159 : C SEGACT MPOVOL*MOD 160 : C SEGACT MPNORM*MOD 161 : C 161 : CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC) 162 : C SEGACT MPORES*MOD 163 : C 164 : C 164 : CALL LICHT(IRN,MPRN,TYPE,IGEOMC) 164 : CALL LICHT(IGN,MPGN,TYPE,IGEOMC) 164 : CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC) 164 : CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC) 164 : CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM) 165 : C 166 : C SEGACT MPRN*MOD 167 : C SEGACT MPGN*MOD 168 : C SEGACT MPRETN*MOD 169 : C SEGACT MPGAMN*MOD 170 : C SEGACT MPLIM*MOD 171 : C 172 : IF(IGEOMF .NE. MELEMF)THEN 173 : WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso' 174 : WRITE(IOIMP,*) 'Probleme de SPG' 175 : C 21 2 176 : C Données incompatibles 176 : CALL ERREUR(21) 177 : GOTO 9999 178 : ENDIF 179 : IF(IGEOMC .NE. MELEMC)THEN 180 : WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso' 181 : WRITE(IOIMP,*) 'Probleme de SPG' 182 : C 21 2 183 : C Données incompatibles 183 : CALL ERREUR(21) 184 : GOTO 9999 185 : ENDIF 186 : C 187 : C 188 : C**** BOUCLE SUR FACEL pour le calcul du FLUX 189 : C 190 : DO NLCF = 1, NFAC 191 : C 192 : C******* NLCF = numero local du centre de facel 193 : C NGCF = numero global du centre de facel 194 : C NLCF1 = numero local du centre de face 195 : C NGCEG = numero global du centre ELT "gauche" 196 : C NLCEG = numero local du centre ELT "gauche" 197 : C NGCED = numero global du centre ELT "droite" 198 : C NLCED = numero local du centre ELT "droite" 199 : C 200 : NGCEG = IPT2.NUM(1,NLCF) 201 : NGCED = IPT2.NUM(3,NLCF) 202 : NGCF = IPT2.NUM(2,NLCF) 203 : NLCF1 = MLENT2.LECT(NGCF) 204 : NLCEG = MLENT1.LECT(NGCEG) 205 : NLCED = MLENT1.LECT(NGCED) 206 : C 207 : C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris. 208 : C 209 : IF(NLCF .NE. NLCF1)THEN 210 : WRITE(IOIMP,*) 'Anomalie dedans kfmm1.eso' 211 : WRITE(IOIMP,*) 'Probleme de SPG' 212 : C 21 2 213 : C Données incompatibles 213 : CALL ERREUR(21) 214 : GOTO 9999 215 : ENDIF 216 : C 217 : CNX = MPNORM.VPOCHA(NLCF,1) 218 : CNY = MPNORM.VPOCHA(NLCF,2) 219 : CTX = -1.0D0 * CNY 220 : CTY = CNX 221 : C 222 : C 223 : C******* Recuperation des Etats "gauche" et "droite" 224 : C 225 : ROG = MPRN.VPOCHA(NLCEG,1) 226 : RUXG = MPGN.VPOCHA(NLCEG,1) 227 : RUYG = MPGN.VPOCHA(NLCEG,2) 228 : UXG = RUXG / ROG 229 : UYG = RUYG / ROG 230 : UNG = (UXG * CNX) + (UYG * CNY) 231 : UTG = (UXG * CTX) + (UYG * CTY) 232 : RETG = MPRETN.VPOCHA(NLCEG,1) 233 : GAMG = MPGAMN.VPOCHA(NLCEG,1) 234 : REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG)) 235 : REC = REC / ROG 236 : PG = (GAMG - 1.0D0) * (RETG - REC) 237 : VOLG = MPOVOL.VPOCHA(NLCEG,1) 238 : C 239 : IF(NLCEG .NE. NLCED)THEN 240 : ROD = MPRN.VPOCHA(NLCED,1) 241 : RUXD = MPGN.VPOCHA(NLCED,1) 242 : RUYD = MPGN.VPOCHA(NLCED,2) 243 : UXD = RUXD / ROD 244 : UYD = RUYD / ROD 245 : UND = (UXD * CNX) + (UYD * CNY) 246 : RETD = MPRETN.VPOCHA(NLCED,1) 247 : GAMD = MPGAMN.VPOCHA(NLCED,1) 248 : REC = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD)) 249 : REC = REC / ROD 250 : PD = (GAMD - 1.0D0) * (RETD - REC) 251 : VOLD = MPOVOL.VPOCHA(NLCED,1) 252 : ELSE 253 : C Murs au condition aux limite 254 : NLLIM=MLELIM.LECT(NGCF) 255 : IF(NLLIM .EQ.0)THEN 256 : C Mur 257 : ROD = ROG 258 : UND = -1.0D0 * UNG 259 : UTD = UTG 260 : UXD = UND * CNX + UTD * CTX 261 : UYD = UND * CNY + UTD * CTY 262 : RUXD = UXD * ROD 263 : RUYD = UYD * ROD 264 : GAMD = GAMG 265 : RETD = RETG 266 : PD = PG 267 : VOLD = VOLG 268 : ELSE 269 : C Boundary 270 : ROD = MPLIM.VPOCHA(NLLIM,1) 271 : UXD = MPLIM.VPOCHA(NLLIM,2) 272 : UYD = MPLIM.VPOCHA(NLLIM,3) 273 : PD = MPLIM.VPOCHA(NLLIM,4) 274 : GAMD = GAMG 275 : UND = (UXD * CNX) + (UYD * CNY) 276 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD)+ 277 : & (0.5D0*ROD*((UXD*UXD)+(UYD*UYD))) 278 : VOLD = VOLG 279 : ENDIF 280 : ENDIF 281 : C 282 : RHTG=RETG+PG 283 : RHTD=RETD+PD 284 : C 285 : SURF = MPOVSU.VPOCHA(NLCF,1) 286 : SIGMAF = (0.5D0 * (GAMG + GAMD)) * (PG + PD) / (ROG + ROD) 287 : SIGMAF = SIGMAF ** 0.5D0 288 : SIGMAF = (0.5D0 * (ABS(UNG) + ABS(UND))) + SIGMAF 289 : C 290 : FR=(ROG*UNG)+(ROD*UND) 291 : FRUX=(ROG*UNG*UXG)+(ROD*UND*UXD)+((PG+PD)*CNX) 292 : FRUY=(ROG*UNG*UYG)+(ROD*UND*UYD)+((PG+PD)*CNY) 293 : FRET=(UNG*RHTG)+(UND*RHTD) 294 : C 295 : COEF=SURF/(2.0D0*VOLG) 296 : MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) + 297 : & (COEF * (FR - (SIGMAF * ROD))) 298 : MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) + 299 : & (COEF * (FRUX - (SIGMAF * ROD * UXD))) 300 : MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) + 301 : & (COEF * (FRUY - (SIGMAF * ROD * UYD))) 302 : MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) + 303 : & (COEF * (FRET - (SIGMAF * RETD))) 304 : C 305 : IF(NLCEG .NE. NLCED)THEN 306 : C 307 : C********** Domain interieur 308 : C 309 : COEF=SURF/(2.0D0*VOLD) 310 : MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) - 311 : & (COEF * (FR + (SIGMAF * ROG))) 312 : MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) - 313 : & (COEF * (FRUX + (SIGMAF * ROG * UXG))) 314 : MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) - 315 : & (COEF * (FRUY + (SIGMAF * ROG * UYG))) 316 : MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) - 317 : & (COEF * (FRET + (SIGMAF * RETG))) 318 : ENDIF 319 : ENDDO 320 : C 321 : SEGDES MPOVSU 322 : SEGDES MPOVOL 323 : SEGDES MPNORM 324 : C 325 : SEGDES MPORES 326 : C 327 : SEGDES MPRN 328 : SEGDES MPRETN 329 : SEGDES MPGN 330 : SEGDES MPGAMN 331 : C 332 : SEGSUP MLENT1 333 : SEGSUP MLENT2 334 : SEGDES IPT2 335 : C 336 : SEGSUP MLELIM 337 : IF(MPLIM .GT. 0) SEGDES MPLIM 338 : C 339 : C 340 : 9999 CONTINUE 341 : C 342 : RETURN 343 : END 344 : C 345 : 346 : 347 : 348 : 349 : 350 : 351 : 352 : 353 : 354 : 355 : 356 :
© Cast3M 2003 - All rights reserved.
Disclaimer