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