1 : C XLAP13 SOURCE CHAT 05/01/13 04:13:48 5004 2 : SUBROUTINE XLAP13(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 : XLAP13 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,UZG 129 : & ,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG 130 : & ,UXD,UYD,UZD 131 : & ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH 132 : & ,UXF,UYF,UZF,DUXXF,DUXYF,DUXZF,DUYXF,DUYYF,DUYZF 133 : & ,DUZXF,DUZYF,DUZZF 134 : & ,DTXF,DTYF,DTZF 135 : & ,DSTDU,TAUXX,TAUXY,TAUYY,TAUXZ,TAUYZ,TAUZZ 136 : $ ,QX,QY,QZ,XF,YF,ZF 137 : & ,CNX,CNY,CNZ,ORIENT,RO,DIAM,DIAM2,CELL,SURF 138 : C 139 : CHARACTER*8 TYPE 140 : C 141 : C**** Initialisation de 1/DT 142 : C 143 : UNSDT = 0.0D0 144 : C 145 : C**** KRIPAD pour la correspondance global/local de centre 146 : C 146 : CALL KRIPAD(MELEMC,MLCENT) 147 : C 148 : C EN KRIPAD 149 : C SEGACT MELEMC 150 : C SEGACT MLCENT 151 : C 151 : CALL LICHT(IMUC,MPMUC,TYPE,IGEOM) 151 : CALL LICHT(IKAPC,MPKAPC,TYPE,IGEOM) 151 : CALL LICHT(ICVC,MPCVC,TYPE,IGEOM) 151 : CALL LICHT(IROC,MPROC,TYPE,IGEOM) 151 : CALL LICHT(IVITC,MPVITC,TYPE,IGEOM) 151 : CALL LICHT(IGRVC,MPGRVF,TYPE,IGEOM) 151 : CALL LICHT(IGRTC,MPGRTF,TYPE,IGEOM) 151 : CALL LICHT(ISURF,MPSURF,TYPE,IGEOM) 151 : CALL LICHT(INORM,MPNORM,TYPE,IGEOM) 151 : CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM) 151 : CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM) 152 : C 153 : C EN LICHT 154 : C SEGACT*MOD MPMUC 155 : C SEGACT*MOD MPKAPC 156 : C SEGACT*MOD MPCVC 157 : C SEGACT*MOD MPROC 158 : C SEGACT*MOD MPVITC 159 : C SEGACT*MOD MPGRVF 160 : C SEGACT*MOD MPGRTF 161 : C SEGACT*MOD MPSURF 162 : C SEGACT*MOD MPNORM 163 : C SEGACT*MOD MPDIAM 164 : C SEGACT*MOD MPFLUX 165 : C 166 : IF(IVIMP .NE. 0)THEN 166 : CALL LICHT(IVIMP,MPVIMP,TYPE,IGEOM) 167 : C SEGACT*MOD MPVIMP 167 : CALL KRIPAD(IGEOM,MLEVIM) 168 : C SEGACT IGEOM 169 : C SEGACT MLEVIM 170 : MELEME = IGEOM 171 : SEGDES MELEME 172 : ENDIF 173 : IF(ITAUIM .NE. 0)THEN 173 : CALL LICHT(ITAUIM,MPTAUI,TYPE,IGEOM) 174 : C SEGACT*MOD MPTAUI 174 : CALL KRIPAD(IGEOM,MLETAI) 175 : C SEGACT IGEOM 176 : C SEGACT MLETAI 177 : MELEME = IGEOM 178 : SEGDES MELEME 179 : ENDIF 180 : IF(IQIMP .NE. 0)THEN 180 : CALL LICHT(IQIMP,MPQIMP,TYPE,IGEOM) 181 : C SEGACT*MOD MPQIMP 181 : CALL KRIPAD(IGEOM,MLEQIM) 182 : C SEGACT IGEOM 183 : C SEGACT MLEQIM 184 : MELEME = IGEOM 185 : SEGDES MELEME 186 : ENDIF 187 : C 188 : SEGACT MELEFL 189 : SEGACT MELEMF 190 : NFAC = MELEMF.NUM(/2) 191 : C 192 : C**** Boucle sur les faces 193 : C 194 : DO NLCF = 1, NFAC, 1 195 : C 196 : C******* NLCF = numero local du centre de facel 197 : C NGCF = numero global du centre de facel 198 : C NLCF1 = numero local du centre de face 199 : C NGCEG = numero global du centre ELT "gauche" 200 : C NLCEG = numero local du centre ELT "gauche" 201 : C NGCED = numero global du centre ELT "droite" 202 : C NLCED = numero local du centre ELT "droite" 203 : C 204 : NGCF = MELEMF.NUM(1,NLCF) 205 : NGCF1 = MELEFL.NUM(2,NLCF) 206 : IF(NGCF .NE. NGCF1)THEN 207 : MOTERR(1:40)= 'FACEL et FACE = ? ' 207 : CALL ERREUR(5) 208 : GOTO 9999 209 : ENDIF 210 : C 211 : NGCEG = MELEFL.NUM(1,NLCF) 212 : NGCED = MELEFL.NUM(3,NLCF) 213 : NLCEG = MLCENT.LECT(NGCEG) 214 : NLCED = MLCENT.LECT(NGCED) 215 : C 216 : C******* On controlle si sur NGCF on impose de CL 217 : C 218 : C NLFVI = numero local du centre de face sul le maillage des 219 : C "vitesses" "imposées" 220 : C 221 : C NLFTI = numero local du centre de face sul le maillage des 222 : C "tau" "imposés" 223 : C 224 : C NLFQI = numero local du centre de face sul le maillage des 225 : C "q" "imposés" 226 : C 227 : IF(IVIMP .NE. 0)THEN 228 : NLFVI = MLEVIM.LECT(NGCF) 229 : ELSE 230 : NLFVI = 0 231 : ENDIF 232 : C 233 : IF(ITAUIM .NE. 0)THEN 234 : NLFTI = MLETAI.LECT(NGCF) 235 : ELSE 236 : NLFTI = 0 237 : ENDIF 238 : C 239 : IF(IQIMP .NE. 0)THEN 240 : NLFQI = MLEQIM.LECT(NGCF) 241 : ELSE 242 : NLFQI = 0 243 : ENDIF 244 : C 245 : IF(NGCEG .NE. NGCED)THEN 246 : C 247 : C********** Parametres geometriques 248 : C 249 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 250 : XF = MCOORD.XCOOR(ICOORX) 251 : YF = MCOORD.XCOOR(ICOORX+1) 252 : ZF = MCOORD.XCOOR(ICOORX+2) 253 : C 254 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 255 : XG = MCOORD.XCOOR(ICOORX) 256 : YG = MCOORD.XCOOR(ICOORX+1) 257 : ZG = MCOORD.XCOOR(ICOORX+2) 258 : XFMXG = XF - XG 259 : YFMYG = YF - YG 260 : ZFMZG = ZF - ZG 261 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG)) 262 : C 263 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1 264 : XD = MCOORD.XCOOR(ICOORX) 265 : YD = MCOORD.XCOOR(ICOORX+1) 266 : ZD = MCOORD.XCOOR(ICOORX+2) 267 : XFMXD = XF - XD 268 : YFMYD = YF - YD 269 : ZFMZD = ZF - ZD 270 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD)) 271 : C 272 : C********** F=G -> DRG = 0 -> ALPHA = 0 273 : C 274 : ALPHA=DRG/(DRG+DRD) 275 : UMALPH= 1.0D0 - ALPHA 276 : C 277 : C********** Les valeurs à l'interface 278 : C 279 : C DRG=0 -> F=G 280 : C 281 : C 282 : C********** Tenseur de contraintes visqueux 283 : C 284 : DUXXF = MPGRVF.VPOCHA(NLCF,1) 285 : DUXYF = MPGRVF.VPOCHA(NLCF,2) 286 : DUXZF = MPGRVF.VPOCHA(NLCF,3) 287 : DUYXF = MPGRVF.VPOCHA(NLCF,4) 288 : DUYYF = MPGRVF.VPOCHA(NLCF,5) 289 : DUYZF = MPGRVF.VPOCHA(NLCF,6) 290 : DUZXF = MPGRVF.VPOCHA(NLCF,7) 291 : DUZYF = MPGRVF.VPOCHA(NLCF,8) 292 : DUZZF = MPGRVF.VPOCHA(NLCF,9) 293 : C 294 : IF (NLFTI .GT. 0) THEN 295 : TAUXX = MPTAUI.VPOCHA(NLFTI,1) 296 : TAUYY = MPTAUI.VPOCHA(NLFTI,2) 297 : TAUZZ = MPTAUI.VPOCHA(NLFTI,3) 298 : TAUXY = MPTAUI.VPOCHA(NLFTI,4) 299 : TAUXZ = MPTAUI.VPOCHA(NLFTI,5) 300 : TAUYZ = MPTAUI.VPOCHA(NLFTI,6) 301 : ELSE 302 : MU=UMALPH*MPMUC.VPOCHA(NLCEG,1) + 303 : & ALPHA*MPMUC.VPOCHA(NLCED,1) 304 : DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF + DUZZF) 305 : TAUXX = MU * (2.0D0 * DUXXF - DSTDU) 306 : TAUYY = MU * (2.0D0 * DUYYF - DSTDU) 307 : TAUZZ = MU * (2.0D0 * DUZZF - DSTDU) 308 : TAUXY = MU * (DUXYF + DUYXF) 309 : TAUXZ = MU * (DUXZF + DUZXF) 310 : TAUYZ = MU * (DUZYF + DUYZF) 311 : ENDIF 312 : C 313 : C********** Vitesse 314 : C 315 : IF( NLFVI .GT. 0) THEN 316 : UXF = MPVIMP.VPOCHA(NLFVI,1) 317 : UYF = MPVIMP.VPOCHA(NLFVI,2) 318 : UZF = MPVIMP.VPOCHA(NLFVI,3) 319 : ELSE 320 : C************* Gauche 321 : UXG = MPVITC.VPOCHA(NLCEG,1) 322 : UYG = MPVITC.VPOCHA(NLCEG,2) 323 : UZG = MPVITC.VPOCHA(NLCEG,3) 324 : C************* Droite 325 : UXD = MPVITC.VPOCHA(NLCED,1) 326 : UYD = MPVITC.VPOCHA(NLCED,2) 327 : UZD = MPVITC.VPOCHA(NLCED,3) 328 : C************* Face 329 : UXF = UMALPH * UXG + ALPHA * UXD 330 : UYF = UMALPH * UYG + ALPHA * UYD 331 : UZF = UMALPH * UZG + ALPHA * UZD 332 : C************* Correction de la vitesse lineaire exacte 333 : UXF = UXF + 334 : & (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) + 335 : & (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) + 336 : & (DUXZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA))) 337 : UYF = UYF + 338 : & (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) + 339 : & (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) + 340 : & (DUYZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA))) 341 : UZF = UZF + 342 : & (DUZXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) + 343 : & (DUZYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) + 344 : & (DUZZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA))) 345 : ENDIF 346 : C 347 : C********** Flux de chaleur 348 : C 349 : IF(NLFQI .GT. 0)THEN 350 : QX = MPQIMP.VPOCHA(NLFQI,1) 351 : QY = MPQIMP.VPOCHA(NLFQI,2) 352 : QZ = MPQIMP.VPOCHA(NLFQI,3) 353 : ELSE 354 : C************* Gauche 355 : DTXF = MPGRTF.VPOCHA(NLCF,1) 356 : DTYF = MPGRTF.VPOCHA(NLCF,2) 357 : DTZF = MPGRTF.VPOCHA(NLCF,3) 358 : C 359 : KAPPA=UMALPH*MPKAPC.VPOCHA(NLCEG,1) + 360 : & ALPHA*MPKAPC.VPOCHA(NLCED,1) 361 : QX = -1.0D0 * KAPPA * DTXF 362 : QY = -1.0D0 * KAPPA * DTYF 363 : QZ = -1.0D0 * KAPPA * DTZF 364 : C 365 : ENDIF 366 : ELSE 367 : C 368 : C********** MURS 369 : C 370 : C Etat a gauche = Etat droite 371 : C 372 : ALPHA=0.0D0 373 : UMALPH=1.0D0 374 : C 375 : C********** Parametres geometriques 376 : C 377 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 378 : XF = MCOORD.XCOOR(ICOORX) 379 : YF = MCOORD.XCOOR(ICOORX+1) 380 : ZF = MCOORD.XCOOR(ICOORX+2) 381 : C 382 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 383 : XG = MCOORD.XCOOR(ICOORX) 384 : YG = MCOORD.XCOOR(ICOORX+1) 385 : ZG = MCOORD.XCOOR(ICOORX+2) 386 : XFMXG = XF - XG 387 : YFMYG = YF - YG 388 : ZFMZG = ZF - ZG 389 : C 390 : C********** Tenseur de contraintes visqueux 391 : C 392 : DUXXF = MPGRVF.VPOCHA(NLCF,1) 393 : DUXYF = MPGRVF.VPOCHA(NLCF,2) 394 : DUXZF = MPGRVF.VPOCHA(NLCF,3) 395 : DUYXF = MPGRVF.VPOCHA(NLCF,4) 396 : DUYYF = MPGRVF.VPOCHA(NLCF,5) 397 : DUYZF = MPGRVF.VPOCHA(NLCF,6) 398 : DUZXF = MPGRVF.VPOCHA(NLCF,7) 399 : DUZYF = MPGRVF.VPOCHA(NLCF,8) 400 : DUZZF = MPGRVF.VPOCHA(NLCF,9) 401 : C 402 : IF (NLFTI .GT. 0) THEN 403 : TAUXX = MPTAUI.VPOCHA(NLFTI,1) 404 : TAUYY = MPTAUI.VPOCHA(NLFTI,2) 405 : TAUZZ = MPTAUI.VPOCHA(NLFTI,3) 406 : TAUXY = MPTAUI.VPOCHA(NLFTI,4) 407 : TAUXZ = MPTAUI.VPOCHA(NLFTI,5) 408 : TAUYZ = MPTAUI.VPOCHA(NLFTI,6) 409 : ELSE 410 : MU=UMALPH*MPMUC.VPOCHA(NLCEG,1) + 411 : & ALPHA*MPMUC.VPOCHA(NLCED,1) 412 : DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF) 413 : TAUXX = MU * (2.0D0 * DUXXF - DSTDU) 414 : TAUYY = MU * (2.0D0 * DUYYF - DSTDU) 415 : TAUZZ = MU * (2.0D0 * DUZZF - DSTDU) 416 : TAUXY = MU * (DUXYF + DUYXF) 417 : TAUXZ = MU * (DUXZF + DUZXF) 418 : TAUYZ = MU * (DUZYF + DUYZF) 419 : ENDIF 420 : C 421 : C********** Vitesse 422 : C 423 : IF( NLFVI .GT. 0) THEN 424 : UXF = MPVIMP.VPOCHA(NLFVI,1) 425 : UYF = MPVIMP.VPOCHA(NLFVI,2) 426 : UZF = MPVIMP.VPOCHA(NLFVI,3) 427 : ELSE 428 : UXF = MPVITC.VPOCHA(NLCEG,1) 429 : UYF = MPVITC.VPOCHA(NLCEG,2) 430 : UZF = MPVITC.VPOCHA(NLCEG,3) 431 : C************* Correction de la vitesse lineaire exacte 432 : UXF = UXF + 433 : & (DUXXF * XFMXG ) + 434 : & (DUXYF * YFMYG ) + 435 : & (DUXZF * ZFMZG ) 436 : UYF = UYF + 437 : & (DUYXF * XFMXG ) + 438 : & (DUYYF * YFMYG ) + 439 : & (DUYZF * ZFMZG ) 440 : UZF = UZF + 441 : & (DUZXF * XFMXG ) + 442 : & (DUZYF * YFMYG ) + 443 : & (DUZZF * ZFMZG ) 444 : ENDIF 445 : C 446 : C********** Flux de chaleur 447 : C 448 : IF(NLFQI .GT. 0)THEN 449 : QX = MPQIMP.VPOCHA(NLFQI,1) 450 : QY = MPQIMP.VPOCHA(NLFQI,2) 451 : QZ = MPQIMP.VPOCHA(NLFQI,3) 452 : ELSE 453 : C************* Gauche 454 : DTXF = MPGRTF.VPOCHA(NLCF,1) 455 : DTYF = MPGRTF.VPOCHA(NLCF,2) 456 : DTZF = MPGRTF.VPOCHA(NLCF,3) 457 : C 458 : KAPPA=UMALPH*MPKAPC.VPOCHA(NLCEG,1) + 459 : & ALPHA*MPKAPC.VPOCHA(NLCED,1) 460 : QX = -1.0D0 * KAPPA * DTXF 461 : QY = -1.0D0 * KAPPA * DTYF 462 : QZ = -1.0D0 * KAPPA * DTZF 463 : C 464 : ENDIF 465 : ENDIF 466 : C 467 : C******* On calcule le sign du pruduit scalare 468 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 469 : C 470 : CNX = MPNORM.VPOCHA(NLCF,1) 471 : CNY = MPNORM.VPOCHA(NLCF,2) 472 : CNZ = MPNORM.VPOCHA(NLCF,3) 473 : ORIENT = CNX * XFMXG + CNY * YFMYG + CNZ * ZFMZG 474 : ORIENT = SIGN(1.0D0,ORIENT) 475 : IF(ORIENT .NE. 1.0D0)THEN 476 : MOTERR(1:40)= 477 : & 'LAPN , subroutine xlap13.eso. ' 478 : WRITE(IOIMP,*) MOTERR(1:40) 479 : MOTERR(1:40)= 480 : & 'Orientation normales. ' 481 : WRITE(IOIMP,*) MOTERR(1:40) 481 : CALL ERREUR(5) 482 : GOTO 9999 483 : ENDIF 484 : C 485 : C******* Le flux aux interfaces 486 : C 487 : SURF = MPSURF.VPOCHA(NLCF,1) 488 : MPFLUX.VPOCHA(NLCF,1) = ((TAUXX * CNX) + (TAUXY * CNY) + 489 : & (TAUXZ * CNZ)) 490 : & * SURF * (-1.0D0) 491 : MPFLUX.VPOCHA(NLCF,2) = ((TAUXY * CNX) + (TAUYY * CNY) + 492 : & (TAUYZ * CNZ)) 493 : & * SURF * (-1.0D0) 494 : MPFLUX.VPOCHA(NLCF,3) = ((TAUXZ * CNX) + (TAUYZ * CNY) + 495 : & (TAUZZ * CNZ)) 496 : & * SURF * (-1.0D0) 497 : MPFLUX.VPOCHA(NLCF,4) = ( 498 : & ((TAUXX * UXF + TAUXY * UYF + TAUXZ * UZF - QX) * CNX) + 499 : & ((TAUXY * UXF + TAUYY * UYF + TAUYZ * UZF - QY) * CNY) + 500 : & ((TAUXZ * UXF + TAUYZ * UYF + TAUZZ * UZF - QZ) * CNZ)) 501 : & * SURF * (-1.0D0) 502 : C 503 : C****** Le pas de temps 504 : C 505 : CV=UMALPH*MPCVC.VPOCHA(NLCEG,1) + 506 : & ALPHA*MPCVC.VPOCHA(NLCED,1) 507 : RO=UMALPH*MPROC.VPOCHA(NLCEG,1) + 508 : & ALPHA*MPROC.VPOCHA(NLCED,1) 509 : DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) + 510 : & ALPHA*MPDIAM.VPOCHA(NLCED,1) 511 : DIAM2=DIAM*DIAM 512 : CELL = 4.0D0*MU / (DIAM2*RO) 513 : LAMBRO=KAPPA/CV 514 : CELL = MAX(CELL, (4.0D0*LAMBRO/(DIAM2*RO))) 515 : C 516 : IF(CELL .GT. UNSDT)THEN 517 : UNSDT = CELL 518 : ENDIF 519 : C 520 : ENDDO 521 : C 522 : C 523 : DT = 1.0D0 / (UNSDT + XPETIT) 524 : C 525 : SEGDES MELEFL 526 : SEGDES MELEMF 527 : SEGDES MELEMC 528 : SEGDES MPSURF 529 : SEGDES MPNORM 530 : SEGDES MPDIAM 531 : SEGSUP MLCENT 532 : C 533 : SEGDES MPKAPC 534 : SEGDES MPMUC 535 : SEGDES MPCVC 536 : SEGDES MPROC 537 : SEGDES MPVITC 538 : SEGDES MPGRVF 539 : SEGDES MPGRTF 540 : SEGDES MPFLUX 541 : C 542 : IF(IVIMP .NE. 0) THEN 543 : SEGDES MPVIMP 544 : SEGSUP MLEVIM 545 : ENDIF 546 : IF(ITAUIM .NE. 0)THEN 547 : SEGDES MPTAUI 548 : SEGSUP MLETAI 549 : ENDIF 550 : IF(IQIMP .NE. 0)THEN 551 : SEGDES MPQIMP 552 : SEGDES MLEQIM 553 : ENDIF 554 : C 555 : 9999 CONTINUE 556 : RETURN 557 : END 558 : 559 : 560 : 561 : 562 : 563 : 564 : 565 :
© Cast3M 2003 - All rights reserved.
Disclaimer