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