1 : C YLA12T SOURCE PV 08/09/23 21:15:08 6164 2 : C YLAP12 SOURCE GOUNAND 01/09/10 21:16:17 4198 3 : SUBROUTINE YLA12T(IGRTC,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 : YLA12T 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 facesP12 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, 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 : C 137 : C**** KRIPAD pour la correspondance global/local de centre 138 : C 138 : CALL KRIPAD(MELEMC,MLCENT) 139 : C 140 : C EN KRIPAD 141 : C SEGACT MELEMC 142 : C SEGACT MLCENT 143 : C 143 : CALL LICHT(IGRTC,MPGRTF,TYPE,IGEOM) 143 : CALL LICHT(ISURF,MPSURF,TYPE,IGEOM) 143 : CALL LICHT(INORM,MPNORM,TYPE,IGEOM) 143 : CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM) 143 : CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM) 144 : C 145 : C EN LICHT 146 : C SEGACT*MOD MPROC 147 : C SEGACT*MOD MPVITC 148 : C SEGACT*MOD MPGRVF 149 : C SEGACT*MOD MPGRTF 150 : C SEGACT*MOD MPSURF 151 : C SEGACT*MOD MPNORM 152 : C SEGACT*MOD MPDIAM 153 : C SEGACT*MOD MPFLUX 154 : C 155 : IF(IQIMP .NE. 0)THEN 155 : CALL LICHT(IQIMP,MPQIMP,TYPE,IGEOM) 156 : C SEGACT*MOD MPQIMP 156 : CALL KRIPAD(IGEOM,MLEQIM) 157 : C SEGACT IGEOM 158 : C SEGACT MLEQIM 159 : MELEME = IGEOM 160 : SEGDES MELEME 161 : ENDIF 162 : C 163 : SEGACT MELEFL 164 : SEGACT MELEMF 165 : NFAC = MELEMF.NUM(/2) 166 : C 167 : C**** Boucle sur les faces 168 : C 169 : DO NLCF = 1, NFAC, 1 170 : C 171 : C******* NLCF = numero local du centre de facel 172 : C NGCF = numero global du centre de facel 173 : C NLCF1 = numero local du centre de face 174 : C NGCEG = numero global du centre ELT "gauche" 175 : C NLCEG = numero local du centre ELT "gauche" 176 : C NGCED = numero global du centre ELT "droite" 177 : C NLCED = numero local du centre ELT "droite" 178 : C 179 : NGCF = MELEMF.NUM(1,NLCF) 180 : NGCF1 = MELEFL.NUM(2,NLCF) 181 : IF(NGCF .NE. NGCF1)THEN 182 : MOTERR(1:40)= 'FACEL et FACE = ? ' 182 : CALL ERREUR(5) 183 : GOTO 9999 184 : ENDIF 185 : C 186 : NGCEG = MELEFL.NUM(1,NLCF) 187 : NGCED = MELEFL.NUM(3,NLCF) 188 : NLCEG = MLCENT.LECT(NGCEG) 189 : NLCED = MLCENT.LECT(NGCED) 190 : C 191 : C******* On controlle si sur NGCF on impose de CL 192 : C 193 : C NLFVI = numero local du centre de face sul le maillage des 194 : C "vitesses" "imposées" 195 : C 196 : C NLFTI = numero local du centre de face sul le maillage des 197 : C "tau" "imposés" 198 : C 199 : C NLFQI = numero local du centre de face sul le maillage des 200 : C "q" "imposés" 201 : C 202 : IF(IQIMP .NE. 0)THEN 203 : NLFQI = MLEQIM.LECT(NGCF) 204 : ELSE 205 : NLFQI = 0 206 : ENDIF 207 : C 208 : IF(NGCEG .NE. NGCED)THEN 209 : C 210 : C********** Parametres geometriques 211 : C 212 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 213 : XF = MCOORD.XCOOR(ICOORX) 214 : YF = MCOORD.XCOOR(ICOORX+1) 215 : C 216 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 217 : XG = MCOORD.XCOOR(ICOORX) 218 : YG = MCOORD.XCOOR(ICOORX+1) 219 : XFMXG = XF - XG 220 : YFMYG = YF - YG 221 : DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)) 222 : C 223 : ICOORX = ((IDIM + 1) * (NGCED - 1))+1 224 : XD = MCOORD.XCOOR(ICOORX) 225 : YD = MCOORD.XCOOR(ICOORX+1) 226 : XFMXD = XF - XD 227 : YFMYD = YF - YD 228 : DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)) 229 : C 230 : C********** F=G -> DRG = 0 -> ALPHA = 0 231 : ALPHA=DRG/(DRG+DRD) 232 : UMALPH= 1.0D0 - ALPHA 233 : C 234 : C 235 : C********** Les valeurs à l'interface 236 : C 237 : C DRG=0 -> F=G 238 : C 239 : C 240 : C********** Flux de chaleur 241 : C 242 : c IF(NLFQI .GT. 0)THEN 243 : c QX = MPQIMP.VPOCHA(NLFQI,1) 244 : c ELSE 245 : C************* Gauche 246 : DTXF = MPGRTF.VPOCHA(NLCF,1) 247 : C 248 : C CORRECTION 249 : QX = -1.0D0 * DTXF 250 : C 251 : c ENDIF 252 : ELSE 253 : C 254 : C********** MURS 255 : C 256 : C Etat a gauche = Etat droite 257 : C 258 : 259 : ALPHA=0.0D0 260 : UMALPH=1.0D0 261 : C 262 : C********** Parametres geometriques 263 : C 264 : ICOORX = ((IDIM + 1) * (NGCF - 1))+1 265 : XF = MCOORD.XCOOR(ICOORX) 266 : YF = MCOORD.XCOOR(ICOORX+1) 267 : C 268 : ICOORX = ((IDIM + 1) * (NGCEG - 1))+1 269 : XG = MCOORD.XCOOR(ICOORX) 270 : YG = MCOORD.XCOOR(ICOORX+1) 271 : XFMXG = XF - XG 272 : YFMYG = YF - YG 273 : C 274 : C********** Flux de chaleur 275 : C 276 : c IF(NLFQI .GT. 0)THEN 277 : c QX = MPQIMP.VPOCHA(NLFQI,1) 278 : c ELSE 279 : C************* Gauche 280 : DTXF = MPGRTF.VPOCHA(NLCF,1) 281 : c DTYF = MPGRTF.VPOCHA(NLCF,2) 282 : C 283 : C CORRECTION 284 : QX = -1.0D0 * DTXF 285 : C 286 : c ENDIF 287 : ENDIF 288 : C 289 : C******* On calcule le sign du pruduit scalare 290 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 291 : C 292 : c CNX = MPNORM.VPOCHA(NLCF,1) 293 : c CNY = MPNORM.VPOCHA(NLCF,2) 294 : c ORIENT = CNX * XFMXG + CNY * YFMYG 295 : c ORIENT = SIGN(1.0D0,ORIENT) 296 : c IF(ORIENT .NE. 1.0D0)THEN 297 : c MOTERR(1:40)= 298 : c & 'LAPN , subroutine ylap12.eso. ' 299 : c WRITE(IOIMP,*) MOTERR(1:40) 299 : c CALL ERREUR(5) 300 : c GOTO 9999 301 : c ENDIF 302 : C 303 : C******* Le flux aux interfaces 304 : C 305 : SURF = MPSURF.VPOCHA(NLCF,1) 306 : MPFLUX.VPOCHA(NLCF,1) = 307 : & ( QX ) * SURF 308 : C 309 : C****** Le pas de temps 310 : C 311 : c RO=UMALPH + 312 : c & ALPHA 313 : c DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) + 314 : c & ALPHA*MPDIAM.VPOCHA(NLCED,1) 315 : c DIAM2=DIAM*DIAM 316 : c CELL = 4.0D0*MU / (DIAM2*RO) 317 : c CELL = MAX(CELL, (4.0D0/(DIAM2*RO))) 318 : C 319 : c IF(CELL .GT. UNSDT)THEN 320 : c UNSDT = CELL 321 : c ENDIF 322 : C 323 : ENDDO 324 : C 325 : C 326 : c DT = 1.0D0 / (UNSDT + XPETIT) 327 : DT = 1.0D0 328 : C 329 : SEGDES MELEFL 330 : SEGDES MELEMF 331 : SEGDES MELEMC 332 : SEGDES MPSURF 333 : SEGDES MPNORM 334 : SEGDES MPDIAM 335 : SEGSUP MLCENT 336 : C 337 : SEGDES MPGRTF 338 : SEGDES MPFLUX 339 : C 340 : IF(IQIMP .NE. 0)THEN 341 : SEGDES MPQIMP 342 : SEGDES MLEQIM 343 : ENDIF 344 : C 345 : 9999 CONTINUE 346 : RETURN 347 : END 348 : 349 : 350 : 351 : 352 : 353 : 354 : 355 : 356 : 357 :
© Cast3M 2003 - All rights reserved.
Disclaimer