1 : C KSUPG SOURCE CHAT 05/01/13 01:09:19 5004 2 : C KSUPG SOURCE MAGN 97/07/19 21:18:18 2765 3 : SUBROUTINE KSUPG(FN,GR,PG,XYZ,HR,PGSQ,RPG, 4 : & NES,NP,NPG,IAXI,XCOOR, 5 : & WT,WS,HK,PGSK,RPGK,AIRE, 6 : & UMJ,DUMJ,KDEB,KFIN,LRV,NPX,NPGX, 7 : & TN,IPADT,UN,IPADU,NPTD,NEL,ANUK, 8 : & IDCEN,LE, 9 : & AL,AH,AP, 10 : & DTM1,DT,DTT1,DTT2,DIAEL,NUEL) 11 : C-------------------------------------------------------------------- 12 : C Calcul des fonctions tests associées à la formulation Eléments 13 : C Finis Petrov-Galerkin afin de stabiliser les termes de convection. 14 : C On en profite pour calculer des temps caractéristiques. 15 : C-------------------------------------------------------------------- 16 : C Cette subroutine étant écrite en fortran pur, certains arguments 17 : C servent uniquement au dimensionnement des tableaux. 18 : C Afin de doper les calculs, on decoupe la boucle sur les éléments 19 : C par paquets afin de bénéficier à fond de la vectorisation. 20 : C-------------------------------------------------------------------- 21 : C 22 : C-------------------------- 23 : C Paramètre Entrée/Sortie : 24 : C-------------------------- 25 : C 26 : C /S FN : Fonction de forme associé à la transformation géométrique 27 : C /S GR : Gradient de FN dans l'élément de référence 28 : C /S PG : Poids d'intégration associé à chaque point de Gauss 29 : C /S XYZ : Coordonnées des sommets de l'élément 30 : C /S HR : Gradient des fonctions de forme dans l'élément courant 31 : C /S PGSQ : Produit du poids d'intégration par detJ 32 : C /S RPG : Rayon associé à chaque point de Gauss (cas axisymétrique) 33 : C 34 : C E/ NES : Dimension en espace du problème traité (2 en 2D et 3 en 3D) 35 : C E/ NP : Nombre de points support de ddl pour l'élément considéré 36 : C E/ NPG : Nombre de points de Gauss 37 : C E/ IAXI : Précise l'axe de symétrie dans le cas axisymétrique 38 : C (en fait l'axe de symétrie est toujours y -> IAXI=2) 39 : C E/ XCOOR : Coordonnées de l'ensemble des points (voir SMCOORD.INC) 40 : C Petrov-Galerkin choisi pour stabiliser la convection 41 : C Petrov-Galerkin choisi pour stabiliser la convection 42 : C /S HK : Gradient de FN dans l'élément courant 43 : C /S PGSK : Produit du poids d'intégration par detJ 44 : C /S RPGK : Rayon associé à chaque point de Gauss (cas axisymétrique) 45 : C /S AIRE : Volume de chaque élément 46 : C /S UMJ : Vitesse servant à la convection en chaque point de Gauss 47 : C /S DUMJ : Gradient de Vitesse servant en formulation conservative 48 : C E/ KDEB : Indice du début de la fenetre 49 : C E/ KFIN : Indice de fin de la fenetre 50 : C E/ LRV : Dimension de la fenetre 51 : C E/ NPX : Nombre de points maximum par élément 52 : C E/ NPGX : Nombre de points de Gauss maximum 53 : C E/ TN : Inconnue transportée au pas de temps précédent 54 : C E/ IPADT : Correspondance numérotation locale et globale pour les 55 : C points SOMMET : J=LECT(I) : le point numéro I est rangé 56 : C en Jème position pour le CHAMP TN 57 : C E/ UN : Champ de vitesse transportant aux sommets des éléments 58 : C E/ IPADU : Correspondance numérotation locale et globale pour les 59 : C points SOMMET : J=LECT(I) : le point numéro I est rangé 60 : C en Jème position pour le CHAMP UN 61 : C E/ NPTD : Nombre de points sous-tendant le spg associé à UN 62 : C E/ NEL : Nombre de points sous-tendant le spg associé à XXDXDY 63 : C E/ ANUK : Champ contenant le coefficient de diffusion (diffusivité 64 : C pour l'énergie, rapport viscosité sur densité pour qdm) 65 : C E/ IDCEN : Contient le schéma associé au terme convectif 66 : C (1=CENTREE 2=SUPGDC 3=SUPG 4=TVISQUEUX 5=CNG) 67 : C E/ LE : Maillage de la zone considérée 68 : C E/ AL : Première composante de XXDXDY (longueur carac. en x) 69 : C E/ AH : Deuxième composante de XXDXDY (longueur carac. en y) 70 : C E/ AP : Troiième composante de XXDXDY (longueur carac. en z) 71 : C E/ DTM1 : Pas de temps imposé 72 : C E/S DT : Pas de temps 73 : C /S DTT1 : Temps caractéristique associé à la diffusion/convection 74 : C /S DTT2 : Temps caractéristique associé à la diffusion 75 : C /S DIAEL : Longueur carac. de l'élément le plus pénalisant pour dt 76 : C /S NUEL : Numéro de l'élément le plus pénalisant pour dt 77 : C 78 : C-------------------------------------------------------------------- 79 : IMPLICIT INTEGER(I-N) 80 : IMPLICIT REAL*8 (A-H,O-Z) 81 : -INC CCOPTIO 82 : -INC CCREEL 83 : DIMENSION IPADT(*),IPADU(*),LE(NP,NEL) 84 : 85 : DIMENSION XCOOR(*) 86 : DIMENSION UN(NPTD,IDIM),TN(*) 87 : DIMENSION ANUK(LRV) 88 : 89 : DIMENSION FN(NP,NPG),GR(IDIM,NP,NPG),PG(NPG),XYZ(IDIM,NP) 90 : DIMENSION HR(IDIM,NP,NPG),PGSQ(NPG),RPG(NPG) 91 : DIMENSION WT(LRV,NPX,NPGX),WS(LRV,NPX,NPGX),HK(LRV,3,NPX,NPGX) 92 : DIMENSION PGSK(LRV,NPGX),RPGK(LRV,NPGX),AIRE(LRV) 93 : C 94 : PARAMETER (LRV1=64,NPGX1=64) 95 : DIMENSION AL(LRV),AH(LRV),AP(LRV) 96 : DIMENSION UM(LRV1),UMI(LRV1,3),BMI(LRV1),TOI1(LRV1),TOI2(LRV1) 97 : DIMENSION GRAD(LRV1,3,NPGX1),UMJ(LRV,3,NPG),DUMJ(LRV,3,NPG) 98 : DIMENSION UP(LRV1,NPGX1),UPI(LRV1,3,NPGX1) 99 : DIMENSION BM(LRV1,NPGX1),BP(LRV1,NPGX1) 100 : DIMENSION TO1(LRV1,NPGX1),TO2(LRV1,NPGX1) 101 : DIMENSION CXT(LRV1,NPGX1),CYT(LRV1,NPGX1),CZT(LRV1,NPGX1) 102 : DIMENSION CXY(LRV1,NPGX1),CXZ(LRV1,NPGX1),CYZ(LRV1,NPGX1) 103 : DIMENSION AL2(LRV1),AH2(LRV1),AP2(LRV1),XMB(LRV1) 104 : C 105 : CC0 = 1.D0 105 : IF (LRV.NE.LRV1) CALL ARRET(0) 105 : IF (NPG.GT.NPGX1) CALL ARRET(0) 106 : C 107 : C- Pour les éléments de la fenetre KDEB:KFIN initialisations des 108 : C- données associées à l'élément fini utilisé : fonctions de forme, 109 : C- gradients des fonctions de forme, aire de l'élément, poids des 110 : C- points de Gauss et rayon des points de Gauss pour le cas axi. 111 : C 112 : DO 100 K=KDEB,KFIN 113 : KP = K - KDEB + 1 114 : DO 20 I=1,NP 115 : J = LE(I,K) 116 : DO 10 N=1,IDIM 117 : XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N) 118 : 10 CONTINUE 119 : 20 CONTINUE 119 : CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG, 120 : & NES,IDIM,NP,NPG,IAXI,AIRE(KP)) 121 : C 122 : DO 50 L=1,NPG 123 : DO 40 N=1,IDIM 124 : DO 30 I=1,NP 125 : HK(KP,N,I,L) = HR(N,I,L) 126 : 30 CONTINUE 127 : 40 CONTINUE 128 : PGSK(KP,L) = PGSQ(L) 129 : RPGK(KP,L) = RPG(L) 130 : 50 CONTINUE 131 : 100 CONTINUE 132 : C 133 : C- Initialisations des variables 134 : C 135 : DO 120 N=1,IDIM 136 : DO 110 K=KDEB,KFIN 137 : KP = K - KDEB + 1 138 : UMI(KP,N) = XPETIT 139 : 110 CONTINUE 140 : 120 CONTINUE 141 : C 142 : DO 150 L=1,NPG 143 : DO 140 N=1,IDIM 144 : DO 130 K=KDEB,KFIN 145 : KP = K - KDEB + 1 146 : UMJ(KP,N,L) = XPETIT 147 : DUMJ(KP,N,L) = XPETIT 148 : UPI(KP,N,L) = XPETIT 149 : GRAD(KP,N,L) = XPETIT 150 : 130 CONTINUE 151 : 140 CONTINUE 152 : 150 CONTINUE 153 : C 154 : C- Calcul en chaque élément, pour chaque point de Gauss 155 : C- UMJ : Vitesse 156 : C- GRAD : Gradient de l'inconnue scalaire au pas de temps précédent 157 : C- Le gradient n'est à calculer que pour l'option SUPGDC 158 : C 159 : IF (IDCEN.EQ.2) THEN 160 : DO 190 N=1,IDIM 161 : DO 180 L=1,NPG 162 : DO 170 I=1,NP 163 : DO 160 K=KDEB,KFIN 164 : KP = K - KDEB + 1 165 : NF = IPADU(LE(I,K)) 166 : NT = IPADT(LE(I,K)) 167 : UMJ(KP,N,L) = UMJ(KP,N,L) + UN(NF,N)*FN(I,L) 168 : DUMJ(KP,N,L) = DUMJ(KP,N,L) + UN(NF,N)*HK(KP,N,I,L) 169 : GRAD(KP,N,L) = GRAD(KP,N,L) + TN(NT) *HK(KP,N,I,L) 170 : 160 CONTINUE 171 : 170 CONTINUE 172 : 180 CONTINUE 173 : 190 CONTINUE 174 : ELSE 175 : DO 230 N=1,IDIM 176 : DO 220 L=1,NPG 177 : DO 210 I=1,NP 178 : DO 200 K=KDEB,KFIN 179 : KP = K - KDEB + 1 180 : NF = IPADU(LE(I,K)) 181 : UMJ(KP,N,L) = UMJ(KP,N,L) + UN(NF,N)*FN(I,L) 182 : DUMJ(KP,N,L) = DUMJ(KP,N,L) + UN(NF,N)*HK(KP,N,I,L) 183 : 200 CONTINUE 184 : 210 CONTINUE 185 : 220 CONTINUE 186 : 230 CONTINUE 187 : ENDIF 188 : C 189 : C- Calcul pour chaque élément de 190 : C- UMI : Vitesse moyenne 191 : C- pour les options SUPGDC et SUPG 192 : C 193 : IF (IDCEN.EQ.2.OR.IDCEN.EQ.3) THEN 194 : DO 260 L=1,NPG 195 : DO 250 N=1,IDIM 196 : DO 240 K=KDEB,KFIN 197 : KP = K - KDEB + 1 198 : UMI(KP,N) = UMI(KP,N) + UMJ(KP,N,L)*PGSK(KP,L)/AIRE(KP) 199 : 240 CONTINUE 200 : 250 CONTINUE 201 : 260 CONTINUE 202 : ENDIF 203 : C 204 : C 205 : C==================================================================== 206 : C Cas 2D 207 : C==================================================================== 208 : C 209 : C 210 : IF (IDIM.EQ.2) THEN 211 : C 212 : C- Calcul pour chaque élément de 213 : C- UM : Norme de la vitesse moyenne 214 : C- BMI : Inverse d'un temps caractéristique associé à la convection 215 : C- pour l'option SUPGDC et SUPG 216 : C 217 : IF (IDCEN.EQ.2.OR.IDCEN.EQ.3) THEN 218 : DO 1000 K=KDEB,KFIN 219 : KP = K - KDEB + 1 220 : UM(KP) = UMI(KP,1)*UMI(KP,1) + UMI(KP,2)*UMI(KP,2) 221 : UM(KP) = SQRT(UM(KP)) 222 : BMIX = UMI(KP,1) / AL(KP) 223 : BMIY = UMI(KP,2) / AH(KP) 224 : BMI(KP) = BMIX*BMIX + BMIY*BMIY 225 : BMI(KP) = SQRT(BMI(KP)) + XPETIT 226 : 1000 CONTINUE 227 : ENDIF 228 : C 229 : C- Calcul de facteur géométrique associé à chaque élément 230 : C 231 : DO 1010 K=KDEB,KFIN 232 : KP=K-KDEB+1 233 : AL2(KP) = 1.D0 / AL(KP) / AL(KP) 234 : AH2(KP) = 1.D0 / AH(KP) / AH(KP) 235 : XMB(KP) = (AL(KP)+AH(KP)) / 2.D0 236 : 1010 CONTINUE 237 : C 238 : C- Calcul en chaque élément, pour chaque point de Gauss de 239 : C- GRAD : vecteur unitaire porté par le gradient du champ scalaire 240 : C- UP : projection de la vitesse sur la direction donnée par GRAD 241 : C- UPI : vecteur UP*GRAD 242 : C- BM : inverse d'un temps caractéristique basé sur la vitesse 243 : C- BP : idem BP mais en considérant la projection de V sur GRAD 244 : C- pour l'option SUPGDC 245 : C 246 : IF (IDCEN.EQ.2) THEN 247 : DO 1030 L=1,NPG 248 : DO 1020 K=KDEB,KFIN 249 : KP = K - KDEB + 1 250 : AX = GRAD(KP,1,L)*GRAD(KP,1,L) + GRAD(KP,2,L)*GRAD(KP,2,L) 251 : AX = SQRT(AX) + XPETIT 252 : GRAD(KP,1,L) = GRAD(KP,1,L) / AX 253 : GRAD(KP,2,L) = GRAD(KP,2,L) / AX 254 : UP(KP,L) = GRAD(KP,1,L) * UMJ(KP,1,L) 255 : & + GRAD(KP,2,L) * UMJ(KP,2,L) 256 : UPI(KP,1,L) = GRAD(KP,1,L) * UP(KP,L) 257 : UPI(KP,2,L) = GRAD(KP,2,L) * UP(KP,L) 258 : 1020 CONTINUE 259 : 1030 CONTINUE 260 : DO 1050 L=1,NPG 261 : DO 1040 K=KDEB,KFIN 262 : KP = K - KDEB + 1 263 : BMX = UMJ(KP,1,L) / AL(KP) 264 : BMY = UMJ(KP,2,L) / AH(KP) 265 : BPX = UPI(KP,1,L) / AL(KP) 266 : BPY = UPI(KP,2,L) / AH(KP) 267 : BM(KP,L) = BMX*BMX + BMY*BMY 268 : BP(KP,L) = BPX*BPX + BPY*BPY 269 : BM(KP,L) = SQRT(BM(KP,L)) + XPETIT 270 : BP(KP,L) = SQRT(BP(KP,L)) + XPETIT 271 : 1040 CONTINUE 272 : 1050 CONTINUE 273 : ENDIF 274 : C 275 : C----------------------------- 276 : C- DECENTREMENT suivant IDCEN 277 : C----------------------------- 278 : C On calcule dans chaque cas TO1 et TO2 ainsi que le tenseur 279 : C associé à la viscosité numérique afin d'évaluer le pas de 280 : C temps de stabilité pour les schémas explicites. 281 : C 282 : C---------- 283 : C CENTREE : 284 : C---------- 285 : IF (IDCEN.EQ.1) THEN 286 : SI1 = 0.D0 287 : SI2 = 0.D0 288 : DO 1110 L=1,NPG 289 : DO 1100 K=KDEB,KFIN 290 : KP = K - KDEB + 1 291 : TO1(KP,L) = 0.D0 292 : TO2(KP,L) = 0.D0 293 : CXT(KP,L) = 0.D0 294 : CYT(KP,L) = 0.D0 295 : CXY(KP,L) = 0.D0 296 : 1100 CONTINUE 297 : 1110 CONTINUE 298 : C--------- 299 : C SUPGDC : Base théorique dans : A New FE formulation for computational 300 : C fluid dynamics, II Beyond SUPG, HUGHES et al., in Comp.Meth.Appl.Mech. 301 : C Eng., vol 54, pp 341-355 (1986). 302 : C--------- 303 : ELSEIF (IDCEN.EQ.2) THEN 304 : SI1 = 1.D0 305 : SI2 = 1.D0 306 : DO 1130 L=1,NPG 307 : DO 1120 K=KDEB,KFIN 308 : KP = K - KDEB + 1 309 : C 310 : C- Approximation "doublement asymptotique" basée sur la vitesse moyenne 311 : C- HMK : Distance basé sur la vitesse moyenne 312 : C- ALFA : Peclet de maille basé sur la vitesse moyenne 313 : C 314 : HMK = 2.D0 * UM(KP) / BM(KP,L) 315 : ALFA = UM(KP) * HMK / ANUK(KP) / 2.D0 316 : IF (ALFA.GT.3.D0) THEN 317 : AKSI = 1.D0 318 : ELSE 319 : AKSI = ALFA / 3.D0 320 : ENDIF 321 : CCT = AKSI / BM(KP,L) 322 : XMB(KP) = HMK 323 : C 324 : C- Approximation "doublement asymptotique" basée sur la projection de 325 : C- la vitesse sur le gradient du champ scalaire 326 : C- HMK : Distance basée sur la vitesse projetée 327 : C- ALFA : Peclet de maille basé sur la vitesse projetée 328 : C 329 : HPK = 2.D0 * UP(KP,L) / BP(KP,L) 330 : ALFA = UP(KP,L) * HPK / ANUK(KP) / 2.D0 331 : IF (ALFA.GT.3.D0) THEN 332 : AKSI = 1.D0 333 : ELSE 334 : AKSI = ALFA / 3.D0 335 : ENDIF 336 : CCP = AKSI / BP(KP,L) 337 : C 338 : CPT = CCP - CCT 339 : IF (CPT.GE.0.D0) then 340 : CC2 = CPT * CC0 341 : ELSE 342 : CC2 = 0.D0 343 : ENDIF 344 : C 345 : TO1(KP,L) = CCT 346 : TO2(KP,L) = CC2 347 : CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT 348 : & + UPI(KP,1,L)*UPI(KP,1,L)*CC2 349 : CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT 350 : & + UPI(KP,2,L)*UPI(KP,2,L)*CC2 351 : CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT 352 : & + UPI(KP,1,L)*UPI(KP,2,L)*CC2 353 : 354 : 1120 CONTINUE 355 : 1130 CONTINUE 356 : C------- 357 : C SUPG : 358 : C------- 359 : ELSEIF (IDCEN.EQ.3) THEN 360 : SI1 = 1.D0 361 : SI2 = 1.D0 362 : DO 1150 L=1,NPG 363 : DO 1140 K=KDEB,KFIN 364 : KP = K - KDEB + 1 365 : C 366 : C- Approximation "doublement asymptotique" basée sur la vitesse moyenne 367 : C- HMK : Distance basé sur la vitesse moyenne 368 : C- ALFA : Peclet de maille basé sur la vitesse moyenne 369 : C 370 : HMK = 2.D0 * UM(KP) / BMI(KP) 371 : ALFA = UM(KP) * HMK / ANUK(KP) / 2.D0 372 : IF (ALFA.GT.3.D0) THEN 373 : AKSI = 1.D0 374 : ELSE 375 : AKSI = ALFA / 3.D0 376 : ENDIF 377 : CCT = AKSI / BMI(KP) 378 : XMB(KP) = HMK 379 : C 380 : TO1(KP,L) = CCT 381 : TO2(KP,L) = 0.D0 382 : CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT 383 : CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT 384 : CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT 385 : 1140 CONTINUE 386 : 1150 CONTINUE 387 : C------------------- 388 : C Tenseur Visqueux : 389 : C------------------- 390 : ELSEIF (IDCEN.EQ.4) THEN 391 : SI1 =-1.D0 392 : SI2 = 1.D0 393 : DT19 = DTM1 * 0.5D0 394 : DO 1170 L=1,NPG 395 : DO 1160 K=KDEB,KFIN 396 : KP = K - KDEB + 1 397 : TO1(KP,L) = DT19 398 : TO2(KP,L) = 0.D0 399 : CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*DT19 400 : CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*DT19 401 : CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*DT19 402 : 1160 CONTINUE 403 : 1170 CONTINUE 404 : C----------------------------- 405 : C Crank Nicholson généralisé : 406 : C----------------------------- 407 : ELSEIF (IDCEN.EQ.5) THEN 408 : SI1 =-1.D0 409 : SI2 = 1.D0 410 : DT19 = DTM1/6.D0 411 : DO 1190 L=1,NPG 412 : DO 1180 K=KDEB,KFIN 413 : KP = K - KDEB + 1 414 : TO1(KP,L) = DT19 415 : TO2(KP,L) = 0.D0 416 : 1180 CONTINUE 417 : 1190 CONTINUE 418 : ENDIF 419 : C 420 : C--------------------------- 421 : C Pas de temps de stabilité 422 : C--------------------------- 423 : C 424 : C La viscosité utilisée dans l'évaluation des pas de temps de stabilité 425 : C explicite ajoute à la viscosité physique la viscosité numérique : 426 : C DT1 : Correspond à une CFL de 1 et un Peclet de maille de 1 427 : C CFL=1 -> dx=Vdt | 428 : C |-> dt=2D/V2 429 : C Pe=1 -> dx=2D/V | 430 : C DT2 : Correspond au pas de temps de stabilité lié à la diffusion 431 : C FOU=1 -> dt=0.5 dx2/D 432 : C 433 : IF (IDCEN.NE.5) THEN 434 : DO 2010 L=1,NPG 435 : DO 2000 K=KDEB,KFIN 436 : KP = K - KDEB + 1 437 : DT0 = DT 438 : DT1 = 2.D0 / 439 : & ( UMJ(KP,1,L)*UMJ(KP,1,L)/(CXT(KP,L)+ANUK(KP)) 440 : & + UMJ(KP,2,L)*UMJ(KP,2,L)/(CYT(KP,L)+ANUK(KP)) ) 441 : DT2 = 0.5D0 / 442 : & ( (CXT(KP,L)+ANUK(KP))*AL2(KP) 443 : & + (CYT(KP,L)+ANUK(KP))*AH2(KP) ) 444 : IF (DT1.LT.DT) DT=DT1 445 : IF (DT2.LT.DT) DT=DT2 446 : IF (DT.NE.DT0) THEN 447 : DTT1 = DT1 448 : DTT2 = DT2 449 : DIAEL = XMB(KP) 450 : NUEL = K 451 : ENDIF 452 : 2000 CONTINUE 453 : 2010 CONTINUE 454 : ENDIF 455 : C 456 : C--------------------------------------------------- 457 : C Fonction test pour la formulation Petrov-Galerkin 458 : C--------------------------------------------------- 459 : C Ce qui est diffusion numérique en explicite se transforme en 460 : C ajoutant de la viscosité numérique (Tenseurs visqueux et CNG). 461 : C WS : Fonction test pour la partie explicite 462 : C 463 : DO 2050 L=1,NPG 464 : DO 2050 I=1,NP 465 : DO 2050 K=KDEB,KFIN 466 : KP = K - KDEB + 1 467 : WT(KP,I,L) = FN(I,L) + SI1 * 468 : & (TO1(KP,L)*(UMJ(KP,1,L)*HK(KP,1,I,L)+UMJ(KP,2,L)*HK(KP,2,I,L)) 469 : & +TO2(KP,L)*(UPI(KP,1,L)*HK(KP,1,I,L)+UPI(KP,2,L)*HK(KP,2,I,L))) 470 : WS(KP,I,L) = FN(I,L) + SI2 * 471 : & (TO1(KP,L)*(UMJ(KP,1,L)*HK(KP,1,I,L)+UMJ(KP,2,L)*HK(KP,2,I,L)) 472 : & +TO2(KP,L)*(UPI(KP,1,L)*HK(KP,1,I,L)+UPI(KP,2,L)*HK(KP,2,I,L))) 473 : 2050 CONTINUE 474 : C 475 : C 476 : C==================================================================== 477 : C Cas 3D 478 : C==================================================================== 479 : C 480 : C 481 : ELSEIF (IDIM.EQ.3) THEN 482 : IF (IDCEN.EQ.2.OR.IDCEN.EQ.3) THEN 483 : DO 5000 K=KDEB,KFIN 484 : KP = K - KDEB + 1 485 : UM(KP) = UMI(KP,1)*UMI(KP,1) 486 : & + UMI(KP,2)*UMI(KP,2) 487 : & + UMI(KP,3)*UMI(KP,3) 488 : UM(KP) = SQRT(UM(KP)) 489 : BMIX = UMI(KP,1) / AL(KP) 490 : BMIY = UMI(KP,2) / AH(KP) 491 : BMIZ = UMI(KP,3) / AP(KP) 492 : BMI(KP) = BMIX*BMIX + BMIY*BMIY + BMIZ*BMIZ 493 : BMI(KP) = SQRT(BMI(KP)) + XPETIT 494 : 5000 CONTINUE 495 : ENDIF 496 : * 497 : DO 5010 K=KDEB,KFIN 498 : KP=K-KDEB+1 499 : AL2(KP) = 1.D0 / AL(KP) / AL(KP) 500 : AH2(KP) = 1.D0 / AH(KP) / AH(KP) 501 : AP2(KP) = 1.D0 / AP(KP) / AP(KP) 502 : XMB(KP) = (AL(KP)+AH(KP)+AP(KP)) / 3.D0 503 : 5010 CONTINUE 504 : * 505 : IF (IDCEN.EQ.2) THEN 506 : DO 5030 L=1,NPG 507 : DO 5020 K=KDEB,KFIN 508 : KP = K-KDEB+1 509 : AX = GRAD(KP,1,L)*GRAD(KP,1,L) 510 : & + GRAD(KP,2,L)*GRAD(KP,2,L) 511 : & + GRAD(KP,3,L)*GRAD(KP,3,L) 512 : AX = SQRT(AX) + XPETIT 513 : GRAD(KP,1,L) = GRAD(KP,1,L) / AX 514 : GRAD(KP,2,L) = GRAD(KP,2,L) / AX 515 : GRAD(KP,3,L) = GRAD(KP,3,L) / AX 516 : UP(KP,L) = GRAD(KP,1,L) * UMJ(KP,1,L) 517 : & + GRAD(KP,2,L) * UMJ(KP,2,L) 518 : & + GRAD(KP,3,L) * UMJ(KP,3,L) 519 : UPI(KP,1,L) = GRAD(KP,1,L) * UP(KP,L) 520 : UPI(KP,2,L) = GRAD(KP,2,L) * UP(KP,L) 521 : UPI(KP,3,L) = GRAD(KP,3,L) * UP(KP,L) 522 : 5020 CONTINUE 523 : 5030 CONTINUE 524 : DO 5050 L=1,NPG 525 : DO 5040 K=KDEB,KFIN 526 : KP = K-KDEB+1 527 : BMX = UMJ(KP,1,L) / AL(KP) 528 : BMY = UMJ(KP,2,L) / AH(KP) 529 : BMZ = UMJ(KP,3,L) / AP(KP) 530 : BPX = UPI(KP,1,L) / AL(KP) 531 : BPY = UPI(KP,2,L) / AH(KP) 532 : BPZ = UPI(KP,3,L) / AP(KP) 533 : BM(KP,L) = BMX*BMX + BMY*BMY + BMZ*BMZ 534 : BP(KP,L) = BPX*BPX + BPY*BPY + BPZ*BPZ 535 : BM(KP,L) = SQRT(BM(KP,L)) + XPETIT 536 : BP(KP,L) = SQRT(BP(KP,L)) + XPETIT 537 : 5040 CONTINUE 538 : 5050 CONTINUE 539 : ENDIF 540 : C 541 : C----------------------------- 542 : C- DECENTREMENT suivant IDCEN 543 : C----------------------------- 544 : C 545 : C---------- 546 : C CENTREE : 547 : C---------- 548 : IF (IDCEN.EQ.1) THEN 549 : SI1= 0.D0 550 : SI2= 0.D0 551 : DO 5110 L=1,NPG 552 : DO 5100 K=KDEB,KFIN 553 : KP = K - KDEB + 1 554 : TO1(KP,L) = 0.D0 555 : TO2(KP,L) = 0.D0 556 : CXT(KP,L) = 0.D0 557 : CYT(KP,L) = 0.D0 558 : CZT(KP,L) = 0.D0 559 : CXY(KP,L) = 0.D0 560 : CXZ(KP,L) = 0.D0 561 : CYZ(KP,L) = 0.D0 562 : 5100 CONTINUE 563 : 5110 CONTINUE 564 : C--------- 565 : C SUPGDC : 566 : C--------- 567 : ELSEIF (IDCEN.EQ.2) THEN 568 : SI1= 1.D0 569 : SI2= 1.D0 570 : DO 5130 L=1,NPG 571 : DO 5120 K=KDEB,KFIN 572 : KP = K - KDEB + 1 573 : HMK = 2.D0 * UM(KP) / BM(KP,L) 574 : ALFA = UM(KP)*HMK / ANUK(KP) / 2.D0 575 : IF (ALFA.GT.3.D0) THEN 576 : AKSI = 1.D0 577 : ELSE 578 : AKSI = ALFA / 3.D0 579 : ENDIF 580 : CCT = AKSI / BM(KP,L) 581 : XMB(KP) = HMK 582 : * 583 : HPK = 2.D0 * UP(KP,L) / BP(KP,L) 584 : ALFA = UP(KP,L)*HPK / ANUK(KP) / 2.D0 585 : IF (ALFA.GT.3.D0) THEN 586 : AKSI = 1.D0 587 : ELSE 588 : AKSI = ALFA / 3.D0 589 : ENDIF 590 : CCP = AKSI / BP(KP,L) 591 : * 592 : CPT=CCP-CCT 593 : IF (CPT.GE.0.D0) then 594 : CC2 = CPT * CC0 595 : ELSE 596 : CC2 = 0.D0 597 : ENDIF 598 : * 599 : TO1(KP,L) = CCT 600 : TO2(KP,L) = CC2 601 : CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT 602 : & + UPI(KP,1,L)*UPI(KP,1,L)*CC2 603 : CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT 604 : & + UPI(KP,2,L)*UPI(KP,2,L)*CC2 605 : CZT(KP,L) = UMJ(KP,3,L)*UMJ(KP,3,L)*CCT 606 : & + UPI(KP,3,L)*UPI(KP,3,L)*CC2 607 : CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT 608 : & + UPI(KP,1,L)*UPI(KP,2,L)*CC2 609 : CXZ(KP,L) = UMJ(KP,1,L)*UMJ(KP,3,L)*CCT 610 : & + UPI(KP,1,L)*UPI(KP,3,L)*CC2 611 : CYZ(KP,L) = UMJ(KP,2,L)*UMJ(KP,3,L)*CCT 612 : & + UPI(KP,2,L)*UPI(KP,3,L)*CC2 613 : 614 : 5120 CONTINUE 615 : 5130 CONTINUE 616 : C------- 617 : C SUPG : 618 : C------- 619 : ELSEIF (IDCEN.EQ.3) THEN 620 : SI1= 1.D0 621 : SI2= 1.D0 622 : DO 5150 L=1,NPG 623 : DO 5140 K=KDEB,KFIN 624 : KP = K - KDEB + 1 625 : HMK = 2.D0 * UM(KP) / BM(KP,L) 626 : ALFA = UM(KP)*HMK / ANUK(KP) / 2.D0 627 : IF (ALFA.GT.3.D0) THEN 628 : AKSI = 1.D0 629 : ELSE 630 : AKSI = ALFA / 3.D0 631 : ENDIF 632 : CCT = AKSI / BM(KP,L) 633 : XMB(KP) = HMK 634 : TO1(KP,L) = CCT 635 : TO2(KP,L) = 0.D0 636 : CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*CCT 637 : CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*CCT 638 : CZT(KP,L) = UMJ(KP,3,L)*UMJ(KP,3,L)*CCT 639 : CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*CCT 640 : CXZ(KP,L) = UMJ(KP,1,L)*UMJ(KP,3,L)*CCT 641 : CYZ(KP,L) = UMJ(KP,2,L)*UMJ(KP,3,L)*CCT 642 : 5140 CONTINUE 643 : 5150 CONTINUE 644 : C------------------- 645 : C Tenseur Visqueux : 646 : C------------------- 647 : ELSEIF (IDCEN.EQ.4) THEN 648 : SI1 =-1.D0 649 : SI2 = 1.D0 650 : DT19 = DTM1 * 0.5D0 651 : DO 5170 L=1,NPG 652 : DO 5160 K=KDEB,KFIN 653 : KP = K - KDEB + 1 654 : TO1(KP,L) = DT19 655 : TO2(KP,L) = 0.D0 656 : CXT(KP,L) = UMJ(KP,1,L)*UMJ(KP,1,L)*DT19 657 : CYT(KP,L) = UMJ(KP,2,L)*UMJ(KP,2,L)*DT19 658 : CZT(KP,L) = UMJ(KP,3,L)*UMJ(KP,3,L)*DT19 659 : CXY(KP,L) = UMJ(KP,1,L)*UMJ(KP,2,L)*DT19 660 : CXZ(KP,L) = UMJ(KP,1,L)*UMJ(KP,3,L)*DT19 661 : CYZ(KP,L) = UMJ(KP,2,L)*UMJ(KP,3,L)*DT19 662 : 5160 CONTINUE 663 : 5170 CONTINUE 664 : C----------------------------- 665 : C Crank Nicholson généralisé : 666 : C----------------------------- 667 : ELSEIF(IDCEN.EQ.5)THEN 668 : SI1 =-1.D0 669 : SI2 = 1.D0 670 : DT19 = DTM1/6.D0 671 : DO 5190 L=1,NPG 672 : DO 5180 K=KDEB,KFIN 673 : KP = K - KDEB + 1 674 : TO1(KP,L) = DT19 675 : TO2(KP,L) = 0.D0 676 : 5180 CONTINUE 677 : 5190 CONTINUE 678 : ENDIF 679 : C 680 : C--------------------------- 681 : C Pas de temps de stabilité 682 : C--------------------------- 683 : C 684 : IF (IDCEN.NE.5) THEN 685 : DO 5210 L=1,NPG 686 : DO 5200 K=KDEB,KFIN 687 : KP = K - KDEB + 1 688 : DT0 = DT 689 : DT1 = 2.D0 / 690 : & ( UMJ(KP,1,L)*UMJ(KP,1,L)/(CXT(KP,L)+ANUK(KP)) 691 : & + UMJ(KP,2,L)*UMJ(KP,2,L)/(CYT(KP,L)+ANUK(KP)) 692 : & + UMJ(KP,3,L)*UMJ(KP,3,L)/(CZT(KP,L)+ANUK(KP)) ) 693 : 694 : DT2 = 0.5D0 / 695 : & ( (CXT(KP,L)+ANUK(KP))*AL2(KP) 696 : & + (CYT(KP,L)+ANUK(KP))*AH2(KP) 697 : & + (CZT(KP,L)+ANUK(KP))*AP2(KP) ) 698 : 699 : IF (DT1.LT.DT) DT=DT1 700 : IF (DT2.LT.DT) DT=DT2 701 : IF (DT.NE.DT0) THEN 702 : DTT1 = DT1 703 : DTT2 = DT2 704 : DIAEL = XMB(KP) 705 : NUEL = K 706 : ENDIF 707 : 5200 CONTINUE 708 : 5210 CONTINUE 709 : END IF 710 : C 711 : C--------------------------------------------------- 712 : C Fonction test pour la formulation Petrov-Galerkin 713 : C--------------------------------------------------- 714 : C 715 : DO 5220 L=1,NPG 716 : DO 5220 I=1,NP 717 : DO 5220 K=KDEB,KFIN 718 : KP = K - KDEB + 1 719 : WT(KP,I,L)= FN(I,L) + SI1 720 : & * ( TO1(KP,L) * ( UMJ(KP,1,L)*HK(KP,1,I,L) 721 : & + UMJ(KP,2,L)*HK(KP,2,I,L) 722 : & + UMJ(KP,3,L)*HK(KP,3,I,L) ) 723 : & + TO2(KP,L) * ( UPI(KP,1,L)*HK(KP,1,I,L) 724 : & + UPI(KP,2,L)*HK(KP,2,I,L) 725 : & + UPI(KP,3,L)*HK(KP,3,I,L) )) 726 : WS(KP,I,L)= FN(I,L) + SI2 727 : & * ( TO1(KP,L) * ( UMJ(KP,1,L)*HK(KP,1,I,L) 728 : & + UMJ(KP,2,L)*HK(KP,2,I,L) 729 : & + UMJ(KP,3,L)*HK(KP,3,I,L)) 730 : & + TO2(KP,L) * ( UPI(KP,1,L)*HK(KP,1,I,L) 731 : & + UPI(KP,2,L)*HK(KP,2,I,L) 732 : & + UPI(KP,3,L)*HK(KP,3,I,L) )) 733 : 5220 CONTINUE 734 : ENDIF 735 : RETURN 736 : END 737 : 738 : 739 :
© Cast3M 2003 - All rights reserved.
Disclaimer