1 : C ZLAP12 SOURCE CHAT 05/01/13 04:23:33 5004 2 : SUBROUTINE ZLAP12(PROPHY,IROC,ITEMC, 3 : & ITIMP,IRIMP, 4 : & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM, 5 : & ICHFLU,DT) 6 : C*********************************************************************** 7 : C NOM : ZLAP12 8 : C DESCRIPTION : Calcul des flux diffusifs Navier-Stokes2D multi-espèces. 9 : C 10 : C 11 : C LANGAGE : ESOPE 12 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) 13 : C mél : gounand@semt2.smts.cea.fr 14 : C*********************************************************************** 15 : C APPELES (UTIL) : KRIPAD : MELEME -> (num. globale->locale) 16 : C LICHT : Lecture des pointeurs (maillages, valeurs) 17 : C d'un objet de type MCHPOI. 18 : C ERREUR : Gestion des erreurs par GIBI. 19 : C APPELE PAR : ZLAP11 : Chapeau de l'opérateur Gibiane 'LAPN' 20 : C option 'VF'. 21 : C*********************************************************************** 22 : C ENTREES : PROPHY (type PROPHY) : propriétés des espèces 23 : C PROPH2 (type PROPH2) : précond. de PROPHY 24 : C IROC (type MCHPOI) : masse volumique par 25 : C élément. 26 : C ITEMC (type MCHPOI) : température par 27 : C élément. 28 : C ITIMP (type MCHPOI) : CL de 29 : C Dirichlet sur la température. 30 : C IRIMP (type MCHPOI) : CL de 31 : C Dirichlet sur la densité. 32 : C MELEMC (type MELEME) : maillage des centres des 33 : C éléments. 34 : C MELEMF (type MELEME) : maillage des faces des 35 : C éléments. 36 : C MELEFL (type MELEME) : connectivités face-(centre 37 : C gauche, centre droit). 38 : C ISURF (type MCHPOI) : surface des faces. 39 : C INORM (type MCHPOI) : normale aux faces. 40 : C IDIAM (type MCHPOI) : diamètre des éléments. 41 : C ENTREES/SORTIES : - 42 : C SORTIES : ICHFLU (type MCHPOI) : flux diffusif aux 43 : C interfaces. 44 : C DT (type REAL*8) : pas de temps de stabilité 45 : C (Fourier) 46 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 47 : C*********************************************************************** 48 : C VERSION : v1, 25/02/2002, version initiale 49 : C HISTORIQUE : v1, 25/02/2002, création 50 : C HISTORIQUE : 51 : C HISTORIQUE : 52 : C*********************************************************************** 53 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 54 : C en cas de modification de ce sous-programme afin de faciliter 55 : C la maintenance ! 56 : C*********************************************************************** 57 : IMPLICIT INTEGER(I-N) 58 : -INC CCOPTIO 59 : -INC CCREEL 60 : -INC SMCHPOI 61 : -INC SMELEME 62 : -INC SMCOORD 63 : -INC SMLENTI 64 : -INC SMLMOTS 65 : C 66 : POINTEUR IROC.MCHPOI,ITEMC.MCHPOI,ICDIFF.MCHPOI,IGRYKF.MCHPOI 67 : POINTEUR MPROC.MPOVAL,MPTEMC.MPOVAL, 68 : & MPTIMP.MPOVAL,MPRIMP.MPOVAL,MPCDIF.MPOVAL,MPGRYK.MPOVAL, 69 : & MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL, 70 : & MPFLUX.MPOVAL,MPCDIN.MPOVAL 71 : C 72 : POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME 73 : POINTEUR MLCENT.MLENTI,MLETIM.MLENTI,MLERIM.MLENTI 74 : C 75 : INTEGER NESP 76 : SEGMENT PROPHY 77 : CHARACTER*4 NOMESP(NESP+1) 78 : REAL*8 CV(NESP+1) 79 : REAL*8 R(NESP+1) 80 : REAL*8 H0K(NESP+1) 81 : POINTEUR CDIFF(NESP+1).MCHPOI 82 : POINTEUR YK(NESP+1).MCHPOI 83 : POINTEUR GRADYK(NESP+1).MCHPOI 84 : POINTEUR CGRYK(NESP+1).MCHELM 85 : POINTEUR CLYK(NESP+1).MCHPOI 86 : ENDSEGMENT 87 : SEGMENT LIPOVA 88 : POINTEUR MPGRAD(0).MPOVAL 89 : ENDSEGMENT 90 : SEGMENT LIPOV2 91 : POINTEUR MPDIFF(0).MPOVAL 92 : ENDSEGMENT 93 : C 94 : INTEGER ITIMP,IRIMP 95 : & ,ISURF,INORM,IDIAM,ICHFLU 96 : & ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED 97 : & ,NLCEG,NLCED,NLFTI,NLFRI 98 : & , ICOORX, IGEOM 99 : 100 : REAL*8 DT, UNSDT 101 : & ,XG,YG,XFMXG,YFMYG,DRG 102 : & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH 103 : & ,XF,YF 104 : & ,CNX, CNY, ORIENT, DIAM, DIAM2, CELL, SURF 105 : INTEGER IESP 106 : REAL*8 CPK,DK,DYKDX,DYKDY,HKF,CPN,DN,HNF 107 : REAL*8 RHOD,RHOF,RHOG 108 : REAL*8 TEMD,TEMF,TEMG 109 : C 110 : CHARACTER*8 TYPE 111 : C 112 : C**** Initialisation de 1/DT 113 : C 114 : UNSDT = 0.0D0 115 : C 116 : C**** KRIPAD pour la correspondance global/local de centre 117 : C 117 : CALL KRIPAD(MELEMC,MLCENT) 118 : C 119 : C EN KRIPAD 120 : C SEGACT MELEMC 121 : C SEGACT MLCENT 122 : C 122 : CALL LICHT(IROC,MPROC,TYPE,IGEOM) 122 : CALL LICHT(ITEMC,MPTEMC,TYPE,IGEOM) 122 : CALL LICHT(ISURF,MPSURF,TYPE,IGEOM) 122 : CALL LICHT(INORM,MPNORM,TYPE,IGEOM) 122 : CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM) 122 : CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM) 123 : C 124 : C EN LICHT 125 : C SEGACT*MOD MPROC 126 : C SEGACT*MOD MPTEMC 127 : C SEGACT*MOD MPSURF 128 : C SEGACT*MOD MPNORM 129 : C SEGACT*MOD MPDIAM 130 : C SEGACT*MOD MPFLUX 131 : C 132 : SEGINI LIPOVA 133 : SEGINI LIPOV2 134 : SEGACT PROPHY 135 : NESP=PROPHY.CV(/1)-1 136 : DO IESP=1,NESP+1 137 : ICDIFF=PROPHY.CDIFF(IESP) 138 : * In licht SEGACT LIPOV2.MPDIFF(*) 138 : CALL LICHT(ICDIFF,MPCDIF,TYPE,IGEOM) 139 : LIPOV2.MPDIFF(**)=MPCDIF 140 : IGRYKF=PROPHY.GRADYK(IESP) 141 : * In licht SEGACT LIPOVA.MPGRAD(*) 141 : CALL LICHT(IGRYKF,MPGRYK,TYPE,IGEOM) 142 : LIPOVA.MPGRAD(**)=MPGRYK 143 : ENDDO 144 : 145 : IF(ITIMP .NE. 0)THEN 145 : CALL LICHT(ITIMP,MPTIMP,TYPE,IGEOM) 146 : C SEGACT*MOD MPVIMP 146 : CALL KRIPAD(IGEOM,MLETIM) 147 : C SEGACT IGEOM 148 : C SEGACT MLETIM 149 : MELEME = IGEOM 150 : SEGDES MELEME 151 : ENDIF 152 : IF(IRIMP .NE. 0)THEN 152 : CALL LICHT(IRIMP,MPRIMP,TYPE,IGEOM) 153 : C SEGACT*MOD MPVIMP 153 : CALL KRIPAD(IGEOM,MLERIM) 154 : C SEGACT IGEOM 155 : C SEGACT MLERIM 156 : MELEME = IGEOM 157 : SEGDES MELEME 158 : ENDIF 159 : 160 : C 161 : SEGACT MELEFL 162 : SEGACT MELEMF 163 : NFAC = MELEMF.NUM(/2) 164 : C 165 : C**** Boucle sur les faces 166 : C 167 : DO NLCF = 1, NFAC, 1 168 : C 169 : C******* NLCF = numero local du centre de facel 170 : C NGCF = numero global du centre de facel 171 : C NLCF1 = numero local du centre de face 172 : C NGCEG = numero global du centre ELT "gauche" 173 : C NLCEG = numero local du centre ELT "gauche" 174 : C NGCED = numero global du centre ELT "droite" 175 : C NLCED = numero local du centre ELT "droite" 176 : C 177 : NGCF = MELEMF.NUM(1,NLCF) 178 : NGCF1 = MELEFL.NUM(2,NLCF) 179 : IF(NGCF .NE. NGCF1)THEN 180 : MOTERR(1:40)= 'FACEL et FACE = ? ' 180 : CALL ERREUR(5) 181 : GOTO 9999 182 : ENDIF 183 : C 184 : NGCEG = MELEFL.NUM(1,NLCF) 185 : NGCED = MELEFL.NUM(3,NLCF) 186 : NLCEG = MLCENT.LECT(NGCEG) 187 : NLCED = MLCENT.LECT(NGCED) 188 : C 189 : C******* On controlle si sur NGCF on impose de CL 190 : C 191 : C NLFTI = numero local du centre de face sul le maillage des 192 : C "températures" "imposées" 193 : C 194 : IF(ITIMP .NE. 0)THEN 195 : NLFTI = MLETIM.LECT(NGCF) 196 : ELSE 197 : NLFTI = 0 198 : ENDIF 199 : C 200 : IF(IRIMP .NE. 0)THEN 201 : NLFRI = MLERIM.LECT(NGCF) 202 : ELSE 203 : NLFRI = 0 204 : ENDIF 205 : C 206 : IF(NGCEG .NE. NGCED)THEN 207 : C 208 : C********** Parametres geometriques 209 : C 210 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 211 : XF = MCOORD.XCOOR(ICOORX) 212 : YF = MCOORD.XCOOR(ICOORX+1) 213 : C 214 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 215 : XG = MCOORD.XCOOR(ICOORX) 216 : YG = MCOORD.XCOOR(ICOORX+1) 217 : XFMXG = XF - XG 218 : YFMYG = YF - YG 219 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)) 220 : C 221 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1 222 : XD = MCOORD.XCOOR(ICOORX) 223 : YD = MCOORD.XCOOR(ICOORX+1) 224 : XFMXD = XF - XD 225 : YFMYD = YF - YD 226 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)) 227 : C 228 : C********** F=G -> DRG = 0 -> ALPHA = 0 229 : C 230 : ALPHA=DRG/(DRG+DRD) 231 : UMALPH= 1.0D0 - ALPHA 232 : C 233 : C********** Les valeurs à l'interface 234 : C 235 : C DRG=0 -> F=G 236 : C 237 : C 238 : C********** Température 239 : C 240 : IF( NLFTI .GT. 0) THEN 241 : TEMF = MPTIMP.VPOCHA(NLFTI,1) 242 : ELSE 243 : C************* Gauche 244 : TEMG = MPTEMC.VPOCHA(NLCEG,1) 245 : C************* Droite 246 : TEMD = MPTEMC.VPOCHA(NLCED,1) 247 : C************* Face 248 : TEMF = UMALPH * TEMG + ALPHA * TEMD 249 : ENDIF 250 : C 251 : C********** Rho 252 : C 253 : IF( NLFRI .GT. 0) THEN 254 : RHOF = MPRIMP.VPOCHA(NLFRI,1) 255 : ELSE 256 : C************* Gauche 257 : RHOG = MPROC.VPOCHA(NLCEG,1) 258 : C************* Droite 259 : RHOD = MPROC.VPOCHA(NLCED,1) 260 : C************* Face 261 : RHOF = UMALPH * RHOG + ALPHA * RHOD 262 : ENDIF 263 : ELSE 264 : C 265 : C********** MURS 266 : C 267 : C Etat a gauche = Etat droite 268 : C 269 : ALPHA=0.0D0 270 : UMALPH=1.0D0 271 : C 272 : C********** Parametres geometriques 273 : C 274 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 275 : XF = MCOORD.XCOOR(ICOORX) 276 : YF = MCOORD.XCOOR(ICOORX+1) 277 : C 278 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 279 : XG = MCOORD.XCOOR(ICOORX) 280 : YG = MCOORD.XCOOR(ICOORX+1) 281 : XFMXG = XF - XG 282 : YFMYG = YF - YG 283 : C 284 : C********** Température 285 : C 286 : IF( NLFTI .GT. 0) THEN 287 : TEMF = MPTIMP.VPOCHA(NLFTI,1) 288 : ELSE 289 : TEMF = MPTEMC.VPOCHA(NLCEG,1) 290 : ENDIF 291 : C 292 : C********** Rho 293 : C 294 : IF( NLFRI .GT. 0) THEN 295 : RHOF = MPRIMP.VPOCHA(NLFRI,1) 296 : ELSE 297 : RHOF = MPROC.VPOCHA(NLCEG,1) 298 : ENDIF 299 : ENDIF 300 : C 301 : C******* On calcule le sign du pruduit scalare 302 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 303 : C 304 : CNX = MPNORM.VPOCHA(NLCF,1) 305 : CNY = MPNORM.VPOCHA(NLCF,2) 306 : ORIENT = CNX * XFMXG + CNY * YFMYG 307 : ORIENT = SIGN(1.0D0,ORIENT) 308 : IF(ORIENT .NE. 1.0D0)THEN 309 : MOTERR(1:40)= 310 : & 'LAPN , subroutine zlap12.eso. ' 311 : WRITE(IOIMP,*) MOTERR(1:40) 311 : CALL ERREUR(5) 312 : GOTO 9999 313 : ENDIF 314 : C 315 : C******* Les flux aux interfaces 316 : C 317 : SURF = MPSURF.VPOCHA(NLCF,1) 318 : DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) + 319 : & ALPHA*MPDIAM.VPOCHA(NLCED,1) 320 : DIAM2=DIAM*DIAM 321 : MPCDIN=LIPOV2.MPDIFF(NESP+1) 322 : CPN=PROPHY.CV(NESP+1)+PROPHY.R(NESP+1) 323 : HNF=(CPN*TEMF)+PROPHY.H0K(NESP+1) 324 : DO IESP=1,NESP 325 : MPCDIF=LIPOV2.MPDIFF(IESP) 326 : DK=UMALPH*MPCDIF.VPOCHA(NLCEG,1) + 327 : $ ALPHA*MPCDIF.VPOCHA(NLCED,1) 328 : DN=UMALPH*MPCDIN.VPOCHA(NLCEG,1) + 329 : $ ALPHA*MPCDIN.VPOCHA(NLCED,1) 330 : MPGRYK=LIPOVA.MPGRAD(IESP) 331 : DYKDX=MPGRYK.VPOCHA(NLCF,1) 332 : DYKDY=MPGRYK.VPOCHA(NLCF,2) 333 : CPK=PROPHY.CV(IESP)+PROPHY.R(IESP) 334 : HKF=(CPK*TEMF)+PROPHY.H0K(IESP) 335 : * Contribution pour l'espèce IESP 336 : MPFLUX.VPOCHA(NLCF,IDIM+1+IESP)= 337 : $ RHOF*DK*((DYKDX*CNX)+(DYKDY*CNY))*SURF*(-1.D0) 338 : * Contribution pour l'énergie totale 339 : MPFLUX.VPOCHA(NLCF,IDIM+1)= 340 : $ MPFLUX.VPOCHA(NLCF,IDIM+1)+ 341 : $ (RHOF*((HKF*DK)-(HNF*DN)) 342 : $ *((DYKDX*CNX)+(DYKDY*CNY))*SURF*(-1.D0)) 343 : * Le pas de temps 344 : CELL=(4.0D0*DK)/DIAM2 345 : IF(CELL .GT. UNSDT)THEN 346 : UNSDT=CELL 347 : ENDIF 348 : ENDDO 349 : C DO IESP=1,NESP+1 350 : C MPCDIF=LIPOV2.MPDIFF(IESP) 351 : C DK=UMALPH*MPCDIF.VPOCHA(NLCEG,1) + 352 : C $ ALPHA*MPCDIF.VPOCHA(NLCED,1) 353 : C MPGRYK=LIPOVA.MPGRAD(IESP) 354 : C DYKDX=MPGRYK.VPOCHA(NLCF,1) 355 : C DYKDY=MPGRYK.VPOCHA(NLCF,2) 356 : C CPK=PROPHY.CV(IESP)+PROPHY.R(IESP) 357 : C HKF=(CPK*TEMF)+PROPHY.H0K(IESP) 358 : C* Contribution pour l'espèce IESP 359 : C IF ((IESP.GE.1).AND.(IESP.LE.NESP)) THEN 360 : C MPFLUX.VPOCHA(NLCF,IDIM+1+IESP)= 361 : C $ RHOF*DK*((DYKDX*CNX)+(DYKDY*CNY))*SURF*(-1.D0) 362 : C ENDIF 363 : C* Contribution pour l'énergie totale 364 : C MPFLUX.VPOCHA(NLCF,IDIM+1)= 365 : C $ MPFLUX.VPOCHA(NLCF,IDIM+1)+ 366 : C $ (RHOF*HKF*DK*((DYKDX*CNX)+(DYKDY*CNY))*SURF*(-1.D0)) 367 : C* Le pas de temps 368 : C CELL=(4.0D0*DK)/DIAM2 369 : C IF(CELL .GT. UNSDT)THEN 370 : C UNSDT=CELL 371 : C ENDIF 372 : C ENDDO 373 : C 374 : ENDDO 375 : C 376 : C 377 : DT = 1.0D0 / (UNSDT + XPETIT) 378 : C 379 : SEGDES MELEFL 380 : SEGDES MELEMF 381 : SEGDES MELEMC 382 : SEGDES MPSURF 383 : SEGDES MPNORM 384 : SEGDES MPDIAM 385 : SEGSUP MLCENT 386 : C 387 : SEGDES MPROC 388 : SEGDES MPTEMC 389 : SEGDES PROPHY 390 : SEGDES LIPOVA.MPGRAD(*) 391 : SEGSUP LIPOVA 392 : SEGDES LIPOV2.MPDIFF(*) 393 : SEGSUP LIPOV2 394 : SEGDES MPFLUX 395 : C 396 : IF(ITIMP .NE. 0) THEN 397 : SEGDES MPTIMP 398 : SEGSUP MLETIM 399 : ENDIF 400 : C 401 : IF(IRIMP .NE. 0) THEN 402 : SEGDES MPRIMP 403 : SEGSUP MLERIM 404 : ENDIF 405 : C 406 : 9999 CONTINUE 407 : RETURN 408 : END 409 : 410 : 411 : 412 : 413 : 414 : 415 : 416 :
© Cast3M 2003 - All rights reserved.
Disclaimer