1 : C XLAP12 SOURCE CHAT 05/01/13 04:13:45 5004 2 : SUBROUTINE XLAP12(IMUC,IKAPC,ICVC,IROC,IVITC,IGRVC, 3 : & IGRTC,IVIMP,ITAUIM,IQIMP, 4 : & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM, 5 : & ICHFLU,DT) 6 : C 7 : C************************************************************************ 8 : C 9 : C PROJET : CASTEM 2000 10 : C 11 : C NOM : XLAP12 12 : C 13 : C DESCRIPTION : Subroutine appellée par ZLAP11 14 : C 15 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 16 : C 17 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 18 : C 19 : C************************************************************************ 20 : C 21 : C 22 : C ENTRÉES: 23 : C ******* 24 : C 25 : C IMUC : pointeur du CHAMPOINT "CENTRE" 26 : C viscosité dynamique (kg/m^3 * m^2/s dans le SI) 27 : C 28 : C IKAPC : pointeur du CHAMPOINT "CENTRE" 29 : C conductivité thermique (J / (s m K)) 30 : C 31 : C ICVC : pointeur du CHAMPOINT "CENTRE" 32 : C chaleur spécifique à volume constant (J / (kg K)) 33 : C 34 : C IROC : pointeur du CHAMPOINT "CENTRE" densité (kg/m^3) 35 : C 36 : C IVITC : pointeur du CHAMPOINT "CENTRE" vitesse 37 : C 38 : C IGRVC : pointeur du CHAMPOINT "FACE" gradient de vitesse 39 : C 40 : C IGRTC : pointeur du CHAMPOINT "FACE" gradient de 41 : C température 42 : C 43 : C IVIMP : pointeur de CHAMPOINT vitesse imposé (sur des 44 : C points FACE) 45 : C 46 : C ITAUIM : pointeur de CHAMPOINT tenseur de contraintes 47 : C visqueux imposé (sur des points FACE) 48 : C 49 : C IQIMP : pointeur de CHAMPOINT flux de chaleur imposé 50 : C (sur des points FACE) 51 : C 52 : C MELEMC : pointeur du maillage "CENTRE" 53 : C 54 : C MELEMF : pointeur du maillage "FACE" 55 : C 56 : C MELEFL : pointeur du maillage "FACEL" 57 : C 58 : C ISURF : pointeur du CHAMPOINT "FACE" qui contient les 59 : C surfaces de faces 60 : C 61 : C INORM : pointeur du CHAMPOINT "FACE" qui contient les 62 : C normales aux faces 63 : C 64 : C IDIAM : pointeur du CHAMPOINT "CENTRE" qui contient le 65 : C diamètre des elts 66 : C 67 : C 68 : C SORTIES 69 : C ******* 70 : C 71 : C ICHFLU : pointeur du CHAMPOINT "FACE" qui contient les 72 : C flux diffusives aux interface 73 : C 74 : C DT (REAL*8) : pas de temps de stabilité donné par le critère 75 : C de Fourier 76 : C 77 : C*********************************************************************** 78 : C 79 : C**** N.B.: Traitement des conditions aux bords 80 : C 81 : C 'VIMP' : la vitesse imposé n'importe ou! 82 : C 83 : C 'QIMP' : flux de chaleur imposé n'importe ou 84 : C 85 : C 'TAUI' : tenseur de contraintes visqueux imposé n'importe ou 86 : C 87 : C 88 : IMPLICIT INTEGER(I-N) 89 : -INC CCOPTIO 90 : -INC CCREEL 91 : -INC SMCHPOI 92 : -INC SMELEME 93 : -INC SMCOORD 94 : -INC SMLENTI 95 : -INC SMLMOTS 96 : C 97 : POINTEUR MPMUC.MPOVAL,MPKAPC.MPOVAL,MPCVC.MPOVAL, 98 : $ MPROC.MPOVAL, MPVITC.MPOVAL, MPGRVF.MPOVAL, 99 : & MPGRTF.MPOVAL, 100 : & MPVIMP.MPOVAL, MPTAUI.MPOVAL, MPQIMP.MPOVAL, 101 : & MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL, 102 : & MPFLUX.MPOVAL 103 : C 104 : POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME 105 : POINTEUR MLCENT.MLENTI,MLEVIM.MLENTI,MLEQIM.MLENTI,MLETAI.MLENTI 106 : C 107 : C**** Variables de COOPTIO 108 : C 109 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 110 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 111 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 112 : C & ,IECHO, IIMPI, IOSPI 113 : C & ,IDIM 114 : CC & ,MCOORD 115 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 116 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 117 : C & ,NORINC,NORVAL,NORIND,NORVAD 118 : C & ,NUCROU, IPSAUV 119 : CC 120 : INTEGER IMUC,IKAPC,ICVC, 121 : $ IROC,IVITC,IGRVC,IGRTC,IVIMP,ITAUIM,IQIMP 122 : & ,ISURF,INORM,IDIAM,ICHFLU 123 : & ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED 124 : & ,NLCEG,NLCED,NLFVI,NLFTI,NLFQI 125 : & , ICOORX, IGEOM 126 : 127 : REAL*8 MU,KAPPA,LAMBRO,CV,DT, UNSDT 128 : & ,UXG,UYG 129 : & ,XG,YG,XFMXG,YFMYG,DRG 130 : & ,UXD,UYD 131 : & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH 132 : & ,UXF,UYF,DUXXF,DUXYF,DUYXF,DUYYF,DTXF,DTYF 133 : & ,DSTDU,TAUXX,TAUXY,TAUYY,QX,QY,XF,YF 134 : & ,CNX, CNY, ORIENT, RO, DIAM, DIAM2, CELL, SURF 135 : C 136 : CHARACTER*8 TYPE 137 : C 138 : C**** Initialisation de 1/DT 139 : C 140 : UNSDT = 0.0D0 141 : C 142 : C**** KRIPAD pour la correspondance global/local de centre 143 : C 143 : CALL KRIPAD(MELEMC,MLCENT) 144 : C 145 : C EN KRIPAD 146 : C SEGACT MELEMC 147 : C SEGACT MLCENT 148 : C 148 : CALL LICHT(IMUC,MPMUC,TYPE,IGEOM) 148 : CALL LICHT(IKAPC,MPKAPC,TYPE,IGEOM) 148 : CALL LICHT(ICVC,MPCVC,TYPE,IGEOM) 148 : CALL LICHT(IROC,MPROC,TYPE,IGEOM) 148 : CALL LICHT(IVITC,MPVITC,TYPE,IGEOM) 148 : CALL LICHT(IGRVC,MPGRVF,TYPE,IGEOM) 148 : CALL LICHT(IGRTC,MPGRTF,TYPE,IGEOM) 148 : CALL LICHT(ISURF,MPSURF,TYPE,IGEOM) 148 : CALL LICHT(INORM,MPNORM,TYPE,IGEOM) 148 : CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM) 148 : CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM) 149 : C 150 : C EN LICHT 151 : C SEGACT*MOD MPMUC 152 : C SEGACT*MOD MPKAPC 153 : C SEGACT*MOD MPCVC 154 : C SEGACT*MOD MPROC 155 : C SEGACT*MOD MPVITC 156 : C SEGACT*MOD MPGRVF 157 : C SEGACT*MOD MPGRTF 158 : C SEGACT*MOD MPSURF 159 : C SEGACT*MOD MPNORM 160 : C SEGACT*MOD MPDIAM 161 : C SEGACT*MOD MPFLUX 162 : C 163 : IF(IVIMP .NE. 0)THEN 163 : CALL LICHT(IVIMP,MPVIMP,TYPE,IGEOM) 164 : C SEGACT*MOD MPVIMP 164 : CALL KRIPAD(IGEOM,MLEVIM) 165 : C SEGACT IGEOM 166 : C SEGACT MLEVIM 167 : MELEME = IGEOM 168 : SEGDES MELEME 169 : ENDIF 170 : IF(ITAUIM .NE. 0)THEN 170 : CALL LICHT(ITAUIM,MPTAUI,TYPE,IGEOM) 171 : C SEGACT*MOD MPTAUI 171 : CALL KRIPAD(IGEOM,MLETAI) 172 : C SEGACT IGEOM 173 : C SEGACT MLETAI 174 : MELEME = IGEOM 175 : SEGDES MELEME 176 : ENDIF 177 : IF(IQIMP .NE. 0)THEN 177 : CALL LICHT(IQIMP,MPQIMP,TYPE,IGEOM) 178 : C SEGACT*MOD MPQIMP 178 : CALL KRIPAD(IGEOM,MLEQIM) 179 : C SEGACT IGEOM 180 : C SEGACT MLEQIM 181 : MELEME = IGEOM 182 : SEGDES MELEME 183 : ENDIF 184 : C 185 : SEGACT MELEFL 186 : SEGACT MELEMF 187 : NFAC = MELEMF.NUM(/2) 188 : C 189 : C**** Boucle sur les faces 190 : C 191 : DO NLCF = 1, NFAC, 1 192 : C 193 : C******* NLCF = numero local du centre de facel 194 : C NGCF = numero global du centre de facel 195 : C NLCF1 = numero local du centre de face 196 : C NGCEG = numero global du centre ELT "gauche" 197 : C NLCEG = numero local du centre ELT "gauche" 198 : C NGCED = numero global du centre ELT "droite" 199 : C NLCED = numero local du centre ELT "droite" 200 : C 201 : NGCF = MELEMF.NUM(1,NLCF) 202 : NGCF1 = MELEFL.NUM(2,NLCF) 203 : IF(NGCF .NE. NGCF1)THEN 204 : MOTERR(1:40)= 'FACEL et FACE = ? ' 204 : CALL ERREUR(5) 205 : GOTO 9999 206 : ENDIF 207 : C 208 : NGCEG = MELEFL.NUM(1,NLCF) 209 : NGCED = MELEFL.NUM(3,NLCF) 210 : NLCEG = MLCENT.LECT(NGCEG) 211 : NLCED = MLCENT.LECT(NGCED) 212 : C 213 : C******* On controlle si sur NGCF on impose de CL 214 : C 215 : C NLFVI = numero local du centre de face sul le maillage des 216 : C "vitesses" "imposées" 217 : C 218 : C NLFTI = numero local du centre de face sul le maillage des 219 : C "tau" "imposés" 220 : C 221 : C NLFQI = numero local du centre de face sul le maillage des 222 : C "q" "imposés" 223 : C 224 : IF(IVIMP .NE. 0)THEN 225 : NLFVI = MLEVIM.LECT(NGCF) 226 : ELSE 227 : NLFVI = 0 228 : ENDIF 229 : C 230 : IF(ITAUIM .NE. 0)THEN 231 : NLFTI = MLETAI.LECT(NGCF) 232 : ELSE 233 : NLFTI = 0 234 : ENDIF 235 : C 236 : IF(IQIMP .NE. 0)THEN 237 : NLFQI = MLEQIM.LECT(NGCF) 238 : ELSE 239 : NLFQI = 0 240 : ENDIF 241 : C 242 : IF(NGCEG .NE. NGCED)THEN 243 : C 244 : C********** Parametres geometriques 245 : C 246 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 247 : XF = MCOORD.XCOOR(ICOORX) 248 : YF = MCOORD.XCOOR(ICOORX+1) 249 : C 250 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 251 : XG = MCOORD.XCOOR(ICOORX) 252 : YG = MCOORD.XCOOR(ICOORX+1) 253 : XFMXG = XF - XG 254 : YFMYG = YF - YG 255 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)) 256 : C 257 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1 258 : XD = MCOORD.XCOOR(ICOORX) 259 : YD = MCOORD.XCOOR(ICOORX+1) 260 : XFMXD = XF - XD 261 : YFMYD = YF - YD 262 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)) 263 : C 264 : C********** F=G -> DRG = 0 -> ALPHA = 0 265 : C 266 : ALPHA=DRG/(DRG+DRD) 267 : UMALPH= 1.0D0 - ALPHA 268 : C 269 : C********** Les valeurs à l'interface 270 : C 271 : C DRG=0 -> F=G 272 : C 273 : C 274 : C********** Tenseur de contraintes visqueux 275 : C 276 : DUXXF = MPGRVF.VPOCHA(NLCF,1) 277 : DUXYF = MPGRVF.VPOCHA(NLCF,2) 278 : DUYXF = MPGRVF.VPOCHA(NLCF,3) 279 : DUYYF = MPGRVF.VPOCHA(NLCF,4) 280 : C 281 : IF (NLFTI .GT. 0) THEN 282 : TAUXX = MPTAUI.VPOCHA(NLFTI,1) 283 : TAUYY = MPTAUI.VPOCHA(NLFTI,2) 284 : TAUXY = MPTAUI.VPOCHA(NLFTI,3) 285 : ELSE 286 : MU=UMALPH*MPMUC.VPOCHA(NLCEG,1) + 287 : & ALPHA*MPMUC.VPOCHA(NLCED,1) 288 : DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF) 289 : TAUXX = MU * (2.0D0 * DUXXF - DSTDU) 290 : TAUXY = MU * (DUXYF + DUYXF) 291 : TAUYY = MU * (2.0D0 * DUYYF - DSTDU) 292 : ENDIF 293 : C 294 : C********** Vitesse 295 : C 296 : IF( NLFVI .GT. 0) THEN 297 : UXF = MPVIMP.VPOCHA(NLFVI,1) 298 : UYF = MPVIMP.VPOCHA(NLFVI,2) 299 : ELSE 300 : C************* Gauche 301 : UXG = MPVITC.VPOCHA(NLCEG,1) 302 : UYG = MPVITC.VPOCHA(NLCEG,2) 303 : C************* Droite 304 : UXD = MPVITC.VPOCHA(NLCED,1) 305 : UYD = MPVITC.VPOCHA(NLCED,2) 306 : C************* Face 307 : UXF = UMALPH * UXG + ALPHA * UXD 308 : UYF = UMALPH * UYG + ALPHA * UYD 309 : C************* Correction de la vitesse lineaire exacte 310 : UXF = UXF + 311 : & (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) + 312 : & (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) 313 : UYF = UYF + 314 : & (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) + 315 : & (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) 316 : ENDIF 317 : C 318 : C********** Flux de chaleur 319 : C 320 : IF(NLFQI .GT. 0)THEN 321 : QX = MPQIMP.VPOCHA(NLFQI,1) 322 : QY = MPQIMP.VPOCHA(NLFQI,2) 323 : ELSE 324 : C************* Gauche 325 : DTXF = MPGRTF.VPOCHA(NLCF,1) 326 : DTYF = MPGRTF.VPOCHA(NLCF,2) 327 : C 328 : KAPPA=UMALPH*MPKAPC.VPOCHA(NLCEG,1) + 329 : & ALPHA*MPKAPC.VPOCHA(NLCED,1) 330 : QX = -1.0D0 * KAPPA * DTXF 331 : QY = -1.0D0 * KAPPA * DTYF 332 : C 333 : ENDIF 334 : ELSE 335 : C 336 : C********** MURS 337 : C 338 : C Etat a gauche = Etat droite 339 : C 340 : ALPHA=0.0D0 341 : UMALPH=1.0D0 342 : C 343 : C********** Parametres geometriques 344 : C 345 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 346 : XF = MCOORD.XCOOR(ICOORX) 347 : YF = MCOORD.XCOOR(ICOORX+1) 348 : C 349 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 350 : XG = MCOORD.XCOOR(ICOORX) 351 : YG = MCOORD.XCOOR(ICOORX+1) 352 : XFMXG = XF - XG 353 : YFMYG = YF - YG 354 : C 355 : C********** Tenseur de contraintes visqueux 356 : C 357 : DUXXF = MPGRVF.VPOCHA(NLCF,1) 358 : DUXYF = MPGRVF.VPOCHA(NLCF,2) 359 : DUYXF = MPGRVF.VPOCHA(NLCF,3) 360 : DUYYF = MPGRVF.VPOCHA(NLCF,4) 361 : C 362 : IF (NLFTI .GT. 0) THEN 363 : TAUXX = MPTAUI.VPOCHA(NLFTI,1) 364 : TAUYY = MPTAUI.VPOCHA(NLFTI,2) 365 : TAUXY = MPTAUI.VPOCHA(NLFTI,3) 366 : ELSE 367 : MU=UMALPH*MPMUC.VPOCHA(NLCEG,1) + 368 : & ALPHA*MPMUC.VPOCHA(NLCED,1) 369 : DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF) 370 : TAUXX = MU * (2.0D0 * DUXXF - DSTDU) 371 : TAUXY = MU * (DUXYF + DUYXF) 372 : TAUYY = MU * (2.0D0 * DUYYF - DSTDU) 373 : ENDIF 374 : C 375 : C********** Vitesse 376 : C 377 : IF( NLFVI .GT. 0) THEN 378 : UXF = MPVIMP.VPOCHA(NLFVI,1) 379 : UYF = MPVIMP.VPOCHA(NLFVI,2) 380 : ELSE 381 : UXF = MPVITC.VPOCHA(NLCEG,1) 382 : UYF = MPVITC.VPOCHA(NLCEG,2) 383 : C************* Correction de la vitesse lineaire exacte 384 : UXF = UXF + 385 : & (DUXXF * XFMXG ) + 386 : & (DUXYF * YFMYG ) 387 : UYF = UYF + 388 : & (DUYXF * XFMXG ) + 389 : & (DUYYF * YFMYG ) 390 : ENDIF 391 : C 392 : C********** Flux de chaleur 393 : C 394 : IF(NLFQI .GT. 0)THEN 395 : QX = MPQIMP.VPOCHA(NLFQI,1) 396 : QY = MPQIMP.VPOCHA(NLFQI,2) 397 : ELSE 398 : C************* Gauche 399 : DTXF = MPGRTF.VPOCHA(NLCF,1) 400 : DTYF = MPGRTF.VPOCHA(NLCF,2) 401 : C 402 : KAPPA=UMALPH*MPKAPC.VPOCHA(NLCEG,1) + 403 : & ALPHA*MPKAPC.VPOCHA(NLCED,1) 404 : QX = -1.0D0 * KAPPA * DTXF 405 : QY = -1.0D0 * KAPPA * DTYF 406 : C 407 : ENDIF 408 : ENDIF 409 : C 410 : C******* On calcule le sign du pruduit scalare 411 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 412 : C 413 : CNX = MPNORM.VPOCHA(NLCF,1) 414 : CNY = MPNORM.VPOCHA(NLCF,2) 415 : ORIENT = CNX * XFMXG + CNY * YFMYG 416 : ORIENT = SIGN(1.0D0,ORIENT) 417 : IF(ORIENT .NE. 1.0D0)THEN 418 : MOTERR(1:40)= 419 : & 'LAPN , subroutine xlap12.eso. ' 420 : WRITE(IOIMP,*) MOTERR(1:40) 420 : CALL ERREUR(5) 421 : GOTO 9999 422 : ENDIF 423 : C 424 : C******* Le flux aux interfaces 425 : C 426 : SURF = MPSURF.VPOCHA(NLCF,1) 427 : MPFLUX.VPOCHA(NLCF,1) = ((TAUXX * CNX) + (TAUXY * CNY)) 428 : & * SURF * (-1.0D0) 429 : MPFLUX.VPOCHA(NLCF,2) = ((TAUXY * CNX) + (TAUYY * CNY)) 430 : & * SURF * (-1.0D0) 431 : MPFLUX.VPOCHA(NLCF,3) = ( 432 : & ((TAUXX * UXF + TAUXY * UYF - QX) * CNX) + 433 : & ((TAUXY * UXF + TAUYY * UYF - QY) * CNY)) 434 : & * SURF * (-1.0D0) 435 : C 436 : C****** Le pas de temps 437 : C 438 : CV=UMALPH*MPCVC.VPOCHA(NLCEG,1) + 439 : & ALPHA*MPCVC.VPOCHA(NLCED,1) 440 : RO=UMALPH*MPROC.VPOCHA(NLCEG,1) + 441 : & ALPHA*MPROC.VPOCHA(NLCED,1) 442 : DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) + 443 : & ALPHA*MPDIAM.VPOCHA(NLCED,1) 444 : DIAM2=DIAM*DIAM 445 : CELL = 4.0D0*MU / (DIAM2*RO) 446 : LAMBRO=KAPPA/CV 447 : CELL = MAX(CELL, (4.0D0*LAMBRO/(DIAM2*RO))) 448 : C 449 : IF(CELL .GT. UNSDT)THEN 450 : UNSDT = CELL 451 : ENDIF 452 : C 453 : ENDDO 454 : C 455 : C 456 : DT = 1.0D0 / (UNSDT + XPETIT) 457 : C 458 : SEGDES MELEFL 459 : SEGDES MELEMF 460 : SEGDES MELEMC 461 : SEGDES MPSURF 462 : SEGDES MPNORM 463 : SEGDES MPDIAM 464 : SEGSUP MLCENT 465 : C 466 : SEGDES MPKAPC 467 : SEGDES MPMUC 468 : SEGDES MPCVC 469 : SEGDES MPROC 470 : SEGDES MPVITC 471 : SEGDES MPGRVF 472 : SEGDES MPGRTF 473 : SEGDES MPFLUX 474 : C 475 : IF(IVIMP .NE. 0) THEN 476 : SEGDES MPVIMP 477 : SEGSUP MLEVIM 478 : ENDIF 479 : IF(ITAUIM .NE. 0)THEN 480 : SEGDES MPTAUI 481 : SEGSUP MLETAI 482 : ENDIF 483 : IF(IQIMP .NE. 0)THEN 484 : SEGDES MPQIMP 485 : SEGDES MLEQIM 486 : ENDIF 487 : C 488 : 9999 CONTINUE 489 : RETURN 490 : END 491 : 492 : 493 : 494 : 495 : 496 : 497 : 498 :
© Cast3M 2003 - All rights reserved.
Disclaimer