1 : C ZLAP13 SOURCE CHAT 05/01/13 04:23:36 5004 2 : SUBROUTINE ZLAP13(PROPHY,IROC,ITEMC, 3 : & ITIMP,IRIMP, 4 : & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM, 5 : & ICHFLU,DT) 6 : C*********************************************************************** 7 : C NOM : ZLAP13 8 : C DESCRIPTION : Calcul des flux diffusifs Navier-Stokes3D 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, 05/03/2002, version initiale 49 : C HISTORIQUE : v1, 05/03/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,ZG,XFMXG,YFMYG,ZFMZG,DRG 102 : & ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH 103 : & ,XF,YF,ZF 104 : & ,CNX,CNY,CNZ,ORIENT,DIAM,DIAM2,CELL,SURF 105 : INTEGER IESP 106 : REAL*8 CPK,DK,DYKDX,DYKDY,DYKDZ,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 : ZF = MCOORD.XCOOR(ICOORX+2) 214 : C 215 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 216 : XG = MCOORD.XCOOR(ICOORX) 217 : YG = MCOORD.XCOOR(ICOORX+1) 218 : ZG = MCOORD.XCOOR(ICOORX+2) 219 : XFMXG = XF - XG 220 : YFMYG = YF - YG 221 : ZFMZG = ZF - ZG 222 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG)) 223 : C 224 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1 225 : XD = MCOORD.XCOOR(ICOORX) 226 : YD = MCOORD.XCOOR(ICOORX+1) 227 : ZD = MCOORD.XCOOR(ICOORX+2) 228 : XFMXD = XF - XD 229 : YFMYD = YF - YD 230 : ZFMZD = ZF - ZD 231 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD)) 232 : C 233 : C********** F=G -> DRG = 0 -> ALPHA = 0 234 : C 235 : ALPHA=DRG/(DRG+DRD) 236 : UMALPH= 1.0D0 - ALPHA 237 : C 238 : C********** Les valeurs à l'interface 239 : C 240 : C DRG=0 -> F=G 241 : C 242 : C 243 : C********** Température 244 : C 245 : IF( NLFTI .GT. 0) THEN 246 : TEMF = MPTIMP.VPOCHA(NLFTI,1) 247 : ELSE 248 : C************* Gauche 249 : TEMG = MPTEMC.VPOCHA(NLCEG,1) 250 : C************* Droite 251 : TEMD = MPTEMC.VPOCHA(NLCED,1) 252 : C************* Face 253 : TEMF = UMALPH * TEMG + ALPHA * TEMD 254 : ENDIF 255 : C 256 : C********** Rho 257 : C 258 : IF( NLFRI .GT. 0) THEN 259 : RHOF = MPRIMP.VPOCHA(NLFRI,1) 260 : ELSE 261 : C************* Gauche 262 : RHOG = MPROC.VPOCHA(NLCEG,1) 263 : C************* Droite 264 : RHOD = MPROC.VPOCHA(NLCED,1) 265 : C************* Face 266 : RHOF = UMALPH * RHOG + ALPHA * RHOD 267 : ENDIF 268 : ELSE 269 : C 270 : C********** MURS 271 : C 272 : C Etat a gauche = Etat droite 273 : C 274 : ALPHA=0.0D0 275 : UMALPH=1.0D0 276 : C 277 : C********** Parametres geometriques 278 : C 279 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 280 : XF = MCOORD.XCOOR(ICOORX) 281 : YF = MCOORD.XCOOR(ICOORX+1) 282 : ZF = MCOORD.XCOOR(ICOORX+2) 283 : C 284 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 285 : XG = MCOORD.XCOOR(ICOORX) 286 : YG = MCOORD.XCOOR(ICOORX+1) 287 : ZG = MCOORD.XCOOR(ICOORX+2) 288 : XFMXG = XF - XG 289 : YFMYG = YF - YG 290 : ZFMZG = ZF - ZG 291 : C 292 : C********** Température 293 : C 294 : IF( NLFTI .GT. 0) THEN 295 : TEMF = MPTIMP.VPOCHA(NLFTI,1) 296 : ELSE 297 : TEMF = MPTEMC.VPOCHA(NLCEG,1) 298 : ENDIF 299 : C 300 : C********** Rho 301 : C 302 : IF( NLFRI .GT. 0) THEN 303 : RHOF = MPRIMP.VPOCHA(NLFRI,1) 304 : ELSE 305 : RHOF = MPROC.VPOCHA(NLCEG,1) 306 : ENDIF 307 : ENDIF 308 : C 309 : C******* On calcule le sign du pruduit scalare 310 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 311 : C 312 : CNX = MPNORM.VPOCHA(NLCF,1) 313 : CNY = MPNORM.VPOCHA(NLCF,2) 314 : CNZ = MPNORM.VPOCHA(NLCF,3) 315 : ORIENT = CNX * XFMXG + CNY * YFMYG + CNZ * ZFMZG 316 : ORIENT = SIGN(1.0D0,ORIENT) 317 : IF(ORIENT .NE. 1.0D0)THEN 318 : MOTERR(1:40)= 319 : & 'LAPN , subroutine zlap13.eso. ' 320 : WRITE(IOIMP,*) MOTERR(1:40) 321 : MOTERR(1:40)= 322 : & 'Orientation normales. ' 323 : WRITE(IOIMP,*) MOTERR(1:40) 323 : CALL ERREUR(5) 324 : GOTO 9999 325 : ENDIF 326 : C 327 : C******* Les flux aux interfaces 328 : C 329 : SURF = MPSURF.VPOCHA(NLCF,1) 330 : DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) + 331 : & ALPHA*MPDIAM.VPOCHA(NLCED,1) 332 : DIAM2=DIAM*DIAM 333 : MPCDIN=LIPOV2.MPDIFF(NESP+1) 334 : CPN=PROPHY.CV(NESP+1)+PROPHY.R(NESP+1) 335 : HNF=(CPN*TEMF)+PROPHY.H0K(NESP+1) 336 : DO IESP=1,NESP 337 : MPCDIF=LIPOV2.MPDIFF(IESP) 338 : DK=UMALPH*MPCDIF.VPOCHA(NLCEG,1) + 339 : $ ALPHA*MPCDIF.VPOCHA(NLCED,1) 340 : DN=UMALPH*MPCDIN.VPOCHA(NLCEG,1) + 341 : $ ALPHA*MPCDIN.VPOCHA(NLCED,1) 342 : MPGRYK=LIPOVA.MPGRAD(IESP) 343 : DYKDX=MPGRYK.VPOCHA(NLCF,1) 344 : DYKDY=MPGRYK.VPOCHA(NLCF,2) 345 : DYKDZ=MPGRYK.VPOCHA(NLCF,3) 346 : CPK=PROPHY.CV(IESP)+PROPHY.R(IESP) 347 : HKF=(CPK*TEMF)+PROPHY.H0K(IESP) 348 : * Contribution pour l'espèce IESP 349 : MPFLUX.VPOCHA(NLCF,IDIM+1+IESP)= 350 : $ RHOF*DK*((DYKDX*CNX)+(DYKDY*CNY)+(DYKDZ*CNZ)) 351 : $ *SURF*(-1.D0) 352 : 353 : * Contribution pour l'énergie totale 354 : MPFLUX.VPOCHA(NLCF,IDIM+1)= 355 : $ MPFLUX.VPOCHA(NLCF,IDIM+1)+ 356 : $ (RHOF*((HKF*DK)-(HNF*DN)) 357 : $ *((DYKDX*CNX)+(DYKDY*CNY)+(DYKDZ*CNZ))*SURF*(-1.D0)) 358 : * Le pas de temps 359 : CELL=(4.0D0*DK)/DIAM2 360 : IF(CELL .GT. UNSDT)THEN 361 : UNSDT=CELL 362 : ENDIF 363 : ENDDO 364 : ENDDO 365 : C 366 : C 367 : DT = 1.0D0 / (UNSDT + XPETIT) 368 : C 369 : SEGDES MELEFL 370 : SEGDES MELEMF 371 : SEGDES MELEMC 372 : SEGDES MPSURF 373 : SEGDES MPNORM 374 : SEGDES MPDIAM 375 : SEGSUP MLCENT 376 : C 377 : SEGDES MPROC 378 : SEGDES MPTEMC 379 : SEGDES PROPHY 380 : SEGDES LIPOVA.MPGRAD(*) 381 : SEGSUP LIPOVA 382 : SEGDES LIPOV2.MPDIFF(*) 383 : SEGSUP LIPOV2 384 : SEGDES MPFLUX 385 : C 386 : IF(ITIMP .NE. 0) THEN 387 : SEGDES MPTIMP 388 : SEGSUP MLETIM 389 : ENDIF 390 : C 391 : IF(IRIMP .NE. 0) THEN 392 : SEGDES MPRIMP 393 : SEGSUP MLERIM 394 : ENDIF 395 : C 396 : 9999 CONTINUE 397 : RETURN 398 : END 399 : 400 : 401 : 402 : 403 : 404 : 405 : 406 :
© Cast3M 2003 - All rights reserved.
Disclaimer