1 : C CFL5 SOURCE KK2000 13/11/08 21:15:17 7860 2 : SUBROUTINE CFL5(ICAS,IPMAIL,MELE,IVAM1,IVAM2,MELV1,MELV2,N2) 3 : *---------------------------------------------------------------------* 4 : * 5 : * calcul du pas de temps CFL 6 : * 7 : * elements poutre,tuyau,barre,cerc,coque2,coque3,dkt 8 : * 9 : * 10 : * entree 11 : * icas : cas à traiter 12 : * = 1 calcul du pas de temps complet ivam1 avec cara 13 : * = 2 calcul du pas de temps lorsque cson est donne ivam2 14 : * = 3 calcul du pas de temps lorsque la taille est donnée ivam1 si cara 15 : * = 4 calcul de la vitesse du son ivam1 donné 16 : * = 5 calcul du parametre de taille ivam1 si cara 17 : * 18 : * ipmail : pointeur vers le maillage a traiter 19 : * mele : numero de l'élément finis dans nomtp 20 : * ivam1 : pointeur vers mptval du cham1 actif 21 : * ivam2 : pointeur vers mptval du cham2 actif 22 : * n2 : nombre de composante en sortie 23 : * 24 : * sortie 25 : * melv1 : melval de la première composante du chamelem resultat 26 : * inactif en sortie 27 : * melv2 : melval de la deuxième composante du chamelem resultat 28 : * inactif en sortie 29 : * 30 : * 31 : *---------------------------------------------------------------------* 32 : IMPLICIT INTEGER(I-N) 33 : IMPLICIT REAL*8(A-H,O-Z) 34 : * 35 : DIMENSION V13(3),V12(3),V23(3) 36 : * 37 : -INC CCOPTIO 38 : -INC CCREEL 39 : -INC CCHAMP 40 : -INC SMCHAML 41 : -INC SMELEME 42 : -INC SMCOORD 43 : * 44 : * 45 : SEGMENT MPTVAL 46 : * ipos pointeur vers la sous zone du mchelm 47 : * nsof ?? 48 : INTEGER IPOS(NS) ,NSOF(NS) 49 : * ival pointeur vers le melval de la composante 50 : * =0 si il n'est pas présente 51 : * ncosou = nbrfac + nbrobl 52 : INTEGER IVAL(NCOSOU) 53 : * continent le type de composante 54 : CHARACTER*16 TYVAL(NCOSOU) 55 : ENDSEGMENT 56 : POINTEUR MPTVA1.MPTVAL,MPTVA2.MPTVAL 57 : * 58 : * 59 : * 60 : MPTVA1 = IVAM1 61 : MPTVA2 = IVAM2 62 : * 63 : * 64 : * branchement en fonction de l'élément fini 65 : * 66 : * 0 5 0 5 0 67 : GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 68 : 2 99,99,99,99,99,99,27,27,29,99,99,99,99,99,99,99,99,99,99,99, 69 : * bar 70 : 4 99,29,99,44,99,46,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 71 : 6 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 72 : 8 99,99,99,29,99,99,99,99,99,99,99,99,27,99,99,99,99,99,99,99, 73 : 1 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 74 : 2 99,99,99,99,99,99,99),MELE 75 : * 76 : * cas de la poutre POUT (mince) et de TIMO (épaisse) 77 : * 78 : 29 CONTINUE 79 : * ================ calcul de la vitesse du son 80 : * le resultat est stocké dans melval avec n1ptel =1 81 : IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN 82 : * recherche des paramètre matériau 83 : * module d'young 84 : MELVA3 = MPTVA1.IVAL(1) 85 : * densite 86 : MELVA4 = MPTVA1.IVAL(3) 87 : SEGACT MELVA3,MELVA4 88 : * 89 : N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2)) 90 : N1PTEL = 1 91 : N2PTEL = 0 92 : N2EL = 0 93 : SEGINI MELVAL 94 : MCSON = MELVAL 95 : * boucle sur les éléments pour calculer la vitesse du son 96 : DO 2903 I=1,N1EL 97 : * on prend les valeurs moyennes sur les éléments 98 : YOU1 = 0.D0 99 : I3 = MIN(I,MELVA3.VELCHE(/2)) 100 : DO 2901 J=1,MELVA3.VELCHE(/1) 101 : YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 102 : 2901 CONTINUE 103 : YOU1 = YOU1 / MELVA3.VELCHE(/1) 104 : * 105 : RO1 = 0.D0 106 : I4 = MIN(I,MELVA4.VELCHE(/2)) 107 : DO 2902 J=1,MELVA4.VELCHE(/1) 108 : RO1 = RO1 + MELVA4.VELCHE(J,I4) 109 : 2902 CONTINUE 110 : RO1 = RO1 / MELVA4.VELCHE(/1) 111 : IF (RO1.EQ.0.D0) THEN 112 : SEGDES MELVA4,MELVA3 113 : SEGSUP MELVAL 113 : CALL ERREUR(855) 114 : RETURN 115 : ENDIF 116 : * 117 : IF (YOU1.EQ.0.D0) THEN 118 : SEGDES MELVA4,MELVA3 119 : SEGSUP MELVAL 119 : CALL ERREUR(856) 120 : RETURN 121 : ENDIF 122 : VELCHE(1,I) = SQRT(YOU1/RO1) 123 : * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 124 : 2903 CONTINUE 125 : SEGDES MELVA4,MELVA3 126 : IF (ICAS.EQ.4) THEN 127 : * cas ou seule la vitesse du son est demandée 128 : MELV2 = 0 129 : MELV1 = MELVAL 130 : SEGDES MELVAL 131 : RETURN 132 : ENDIF 133 : ELSE IF (ICAS.EQ.2) THEN 134 : * recuperation du champ 135 : * SEGACT MPTVA2 136 : MELVA1 = MPTVA2.IVAL(1) 137 : SEGACT MELVA1 138 : MCSON = MELVA1 139 : SEGDES MPTVA2 140 : ENDIF 141 : * ================ paramètre geometrique 142 : * stocké dans un melval mtaille 143 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN 144 : * on verifie que le champ de caractéristiques existe 145 : IF (IVAM1.EQ.0) THEN 146 : MOTERR(1:8)='CARACTER' 146 : CALL ERREUR(145) 147 : RETURN 148 : ENDIF 149 : MELEME = IPMAIL 150 : SEGACT MELEME 151 : N1EL = NUM(/2) 152 : N1PTEL = 1 153 : N2PTEL = 0 154 : N2EL = 0 155 : SEGINI MELVAL 156 : MTAIL1 = MELVAL 156 : CALL CFLTAI(MTAIL1,IPMAIL,MELE) 157 : SEGDES MELEME 158 : * second paramètre geométrique 159 : IOBL = 3 160 : IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0 161 : SEGINI, MELVA2=MELVAL 162 : MTAIL2 = MELVA2 163 : IF (MELE.EQ.29) THEN 164 : * récupération de la section des inerties et de la torsion 165 : * section 166 : MELVA3 = MPTVA1.IVAL(4+IOBL) 167 : * inertie z 168 : MELVA4 = MPTVA1.IVAL(3+IOBL) 169 : * inertie y 170 : MELVA5 = MPTVA1.IVAL(2+IOBL) 171 : * torsion 172 : MELVA6 = MPTVA1.IVAL(1+IOBL) 173 : SEGACT MELVA3,MELVA4,MELVA5,MELVA6 174 : DO 2905 I=1,N1EL 175 : I3 = MIN(I,MELVA3.VELCHE(/2)) 176 : I4 = MIN(I,MELVA4.VELCHE(/2)) 177 : I5 = MIN(I,MELVA5.VELCHE(/2)) 178 : I6 = MIN(I,MELVA6.VELCHE(/2)) 179 : * remaarque les valeurs ne sont pas moyennées sur l'élément 180 : DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) 181 : DUM1 = DUM1 * MELVA3.VELCHE(1,I3) 182 : DUM1=DUM1/MAX(MELVA4.VELCHE(1,I4),MELVA5.VELCHE(1,I5), 183 : & MELVA6.VELCHE(1,I6)) 184 : DUM1=SQRT(DUM1) * 0.10132D0 185 : MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 186 : 2905 CONTINUE 187 : SEGDES MELVA3,MELVA4,MELVA5,MELVA6 188 : ELSEIF (MELE.EQ.84) THEN 189 : * récupération de la section des inerties et de la torsion 190 : * section 191 : MELVA3 = MPTVA1.IVAL(4+IOBL) 192 : * inertie z 193 : MELVA4 = MPTVA1.IVAL(3+IOBL) 194 : * inertie y 195 : MELVA5 = MPTVA1.IVAL(2+IOBL) 196 : * torsion 197 : MELVA6 = MPTVA1.IVAL(1+IOBL) 198 : SEGACT MELVA3,MELVA4,MELVA5,MELVA6 199 : DO 8405 I=1,N1EL 200 : I3 = MIN(I,MELVA3.VELCHE(/2)) 201 : I4 = MIN(I,MELVA4.VELCHE(/2)) 202 : I5 = MIN(I,MELVA5.VELCHE(/2)) 203 : I6 = MIN(I,MELVA6.VELCHE(/2)) 204 : * remarque les valeurs ne sont pas moyennées sur l'élément 205 : DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) 206 : * 0.9 = 2**1.5 / XPI 0.7 coef de securite 207 : DUM2 = MELVA2.VELCHE(1,I)*0.9D0*0.7D0 208 : DUM1 = DUM1 * MELVA3.VELCHE(1,I3) 209 : DUM1=DUM1/MAX(MELVA4.VELCHE(1,I4),MELVA5.VELCHE(1,I5), 210 : & MELVA6.VELCHE(1,I6)) 211 : DUM1=SQRT(DUM1) * 0.10132D0 212 : * correction pour tenir compte du cisaillement Blevins page 175 eq 8.30 213 : MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 + DUM2 214 : 8405 CONTINUE 215 : SEGDES MELVA3,MELVA4,MELVA5,MELVA6 216 : ELSEIF(MELE.EQ.42) THEN 217 : * cas du tuyau 218 : * récupération des caractéristiques 219 : * epai 220 : MELVA3 = MPTVA1.IVAL(1+IOBL) 221 : * rayo 222 : MELVA4 = MPTVA1.IVAL(2+IOBL) 223 : * pour les autres on ne ient pas compte des modification 224 : * qui assouplissent le tuyau donc omega max diminue 225 : SEGACT MELVA3,MELVA4 226 : DO 4205 I=1,N1EL 227 : I3 = MIN(I,MELVA3.VELCHE(/2)) 228 : I4 = MIN(I,MELVA4.VELCHE(/2)) 229 : REXT = MELVA4.VELCHE(1,I4) 230 : RINT = REXT - MELVA3.VELCHE(1,I3) 231 : * remarque les valeurs ne sont pas moyennées sur l'élément 232 : DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) 233 : DUM1 = DUM1 * 2.D0 / (REXT*REXT+RINT*RINT) 234 : DUM1=SQRT(DUM1) * 0.10132D0 235 : MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 236 : 4205 CONTINUE 237 : SEGDES MELVA3,MELVA4 238 : ENDIF 239 : IF (ICAS.EQ.5) THEN 240 : MELV1 = MTAIL1 241 : MELV2 = MTAIL2 242 : SEGDES MELVAL,MELVA2 243 : RETURN 244 : ENDIF 245 : ELSE IF (ICAS.EQ.3) THEN 246 : * recuperation du champ 247 : * SEGACT MPTVA2 248 : MELVA1 = MPTVA2.IVAL(1) 249 : MELVA2 = MPTVA2.IVAL(2) 250 : SEGACT MELVA1,MELVA2 251 : MTAIL1 = MELVA1 252 : MTAIL2 = MELVA2 253 : SEGDES MPTVA2 254 : ENDIF 255 : * ================ pas de temps cfl 256 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN 257 : * recuperation de la vitesse du son 258 : * et du paramètre de taille 259 : MELVA1 = MCSON 260 : MELVA2 = MTAIL1 261 : MELVA3 = MTAIL2 262 : * creation du melval résultat 263 : N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2), 264 : & MELVA3.VELCHE(/2)) 265 : N1PTEL = 1 266 : N2EL = 0 267 : N2PTEL = 0 268 : SEGINI MELVAL 269 : * 270 : DO 2904 I=1,N1EL 271 : I1 = MIN(I,MELVA1.VELCHE(/2)) 272 : I2 = MIN(I,MELVA2.VELCHE(/2)) 273 : I3 = MIN(I,MELVA3.VELCHE(/2)) 274 : DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3)) 275 : VELCHE(1,I)= DUM1 / MELVA1.VELCHE(1,I1) 276 : * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 277 : 2904 CONTINUE 278 : MELV1 = MELVAL 279 : MELV2 = 0 280 : SEGDES MELVAL 281 : IF (ICAS.EQ.1) THEN 282 : SEGSUP MELVA1,MELVA2,MELVA3 283 : ELSE IF (ICAS.EQ.2) THEN 284 : SEGSUP MELVA2,MELVA3 285 : SEGDES MELVA1 286 : ELSE 287 : * icas=3 288 : SEGSUP MELVA1 289 : SEGDES MELVA2,MELVA3 290 : ENDIF 291 : RETURN 292 : ENDIF 293 : * 294 : * 295 : * cas du coq2 296 : * 297 : 44 CONTINUE 298 : * ================ calcul de la vitesse du son 299 : * le resultat est stocké dans melval avec n1ptel =1 300 : IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN 301 : * recherche des paramètre matériau 302 : * module d'young 303 : MELVA3 = MPTVA1.IVAL(1) 304 : * densite 305 : MELVA4 = MPTVA1.IVAL(3) 306 : * nu 307 : MELVA5 = MPTVA1.IVAL(2) 308 : * SEGDES MPTVA1 309 : SEGACT MELVA3,MELVA4,MELVA5 310 : * 311 : N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2), 312 : & MELVA5.VELCHE(/2)) 313 : N1PTEL = 1 314 : N2PTEL = 0 315 : N2EL = 0 316 : SEGINI MELVAL 317 : MCSON = MELVAL 318 : * boucle sur les éléments pour calculer la vitesse du son 319 : DO 4404 I=1,N1EL 320 : * on prend les valeurs moyennes sur les éléments 321 : YOU1 = 0.D0 322 : I3 = MIN(I,MELVA3.VELCHE(/2)) 323 : DO 4401 J=1,MELVA3.VELCHE(/1) 324 : YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 325 : 4401 CONTINUE 326 : YOU1 = YOU1 / MELVA3.VELCHE(/1) 327 : * 328 : RO1 = 0.D0 329 : I4 = MIN(I,MELVA4.VELCHE(/2)) 330 : DO 4402 J=1,MELVA4.VELCHE(/1) 331 : RO1 = RO1 + MELVA4.VELCHE(J,I4) 332 : 4402 CONTINUE 333 : RO1 = RO1 / MELVA4.VELCHE(/1) 334 : * 335 : DNU1 = 0.D0 336 : I5 = MIN(I,MELVA5.VELCHE(/2)) 337 : DO 4403 J=1,MELVA5.VELCHE(/1) 338 : DNU1 = DNU1 + MELVA5.VELCHE(J,I5) 339 : 4403 CONTINUE 340 : DNU1 = DNU1 / MELVA5.VELCHE(/1) 341 : 342 : IF (RO1.EQ.0.D0) THEN 343 : SEGDES MELVA4,MELVA3,MELVA5 344 : SEGSUP MELVAL 344 : CALL ERREUR(855) 345 : RETURN 346 : ENDIF 347 : * 348 : IF (YOU1.EQ.0.D0) THEN 349 : SEGDES MELVA4,MELVA3,MELVA5 350 : SEGSUP MELVAL 350 : CALL ERREUR(856) 351 : RETURN 352 : ENDIF 353 : IF (DNU1.GE.1.D0) THEN 354 : SEGDES MELVA4,MELVA3,MELVA5 355 : SEGSUP MELVAL 355 : CALL ERREUR(36) 356 : RETURN 357 : ENDIF 358 : * 359 : VELCHE(1,I) = SQRT(YOU1/(RO1*(1-DNU1*DNU1))) 360 : * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 361 : 4404 CONTINUE 362 : SEGDES MELVA4,MELVA3,MELVA5 363 : IF (ICAS.EQ.4) THEN 364 : * cas ou seule la vitesse du son est demandée 365 : MELVA2 = 0 366 : MELV1 = MELVAL 367 : SEGDES MELVAL 368 : RETURN 369 : ENDIF 370 : ELSE IF (ICAS.EQ.2) THEN 371 : * recuperation du champ 372 : MELVA1 = MPTVA2.IVAL(1) 373 : SEGACT MELVA1 374 : MCSON = MELVA1 375 : SEGDES MPTVA2 376 : ENDIF 377 : * ================ paramètre geometrique 378 : * stocké dans un melval mtaille 379 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN 380 : MELEME = IPMAIL 381 : SEGACT MELEME 382 : N1EL = NUM(/2) 383 : N1PTEL = 1 384 : N2PTEL = 0 385 : N2EL = 0 386 : SEGINI MELVAL 387 : MTAIL1 = MELVAL 387 : CALL CFLTAI(MTAIL1,IPMAIL,MELE) 388 : SEGDES MELEME 389 : * second paramètre geométrique 390 : IOBL = 3 391 : IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0 392 : SEGINI, MELVA2=MELVAL 393 : MTAIL2 = MELVA2 394 : * récupération de l' epaisseur 395 : MELVA3 = MPTVA1.IVAL(1+IOBL) 396 : SEGACT MELVA3 397 : DO 4406 I=1,N1EL 398 : I3 = MIN(I,MELVA3.VELCHE(/2)) 399 : * remarque les valeurs ne sont pas moyennées sur l'élément 400 : DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) 401 : DUM1 = DUM1 / (3.23D0 * MELVA3.VELCHE(1,I3)) 402 : MELVA2.VELCHE(1,I)=DUM1 403 : 4406 CONTINUE 404 : SEGDES MELVA3 405 : SEGDES MPTVA1 406 : * 407 : IF (ICAS.EQ.5) THEN 408 : MELV1 = MTAIL1 409 : MELV2 = MTAIL2 410 : SEGDES MELVAL 411 : RETURN 412 : ENDIF 413 : ELSE IF (ICAS.EQ.3) THEN 414 : * recuperation du champ 415 : MELVA1 = MPTVA2.IVAL(1) 416 : MELVA2 = MPTVA2.IVAL(2) 417 : SEGACT MELVA1,MELVA2 418 : MTAIL1 = MELVA1 419 : MTAIL2 = MELVA2 420 : SEGDES MPTVA2 421 : ENDIF 422 : * ================ pas de temps cfl 423 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN 424 : * recuperation de la vitesse du son 425 : * et du paramètre de taille 426 : MELVA1 = MCSON 427 : MELVA2 = MTAIL1 428 : MELVA3 = MTAIL2 429 : * creation du melval résultat 430 : N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2)) 431 : N1PTEL = 1 432 : N2EL = 0 433 : N2PTEL = 0 434 : SEGINI MELVAL 435 : * 436 : DO 4405 I=1,N1EL 437 : I1 = MIN(I,MELVA1.VELCHE(/2)) 438 : I2 = MIN(I,MELVA2.VELCHE(/2)) 439 : I3 = MIN(I,MELVA3.VELCHE(/2)) 440 : DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3)) 441 : VELCHE(1,I)=DUM1/MELVA1.VELCHE(1,I1) 442 : * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 443 : 4405 CONTINUE 444 : MELV1 = MELVAL 445 : MELV2 = 0 446 : SEGDES MELVAL 447 : IF (ICAS.EQ.1) THEN 448 : SEGSUP MELVA1,MELVA2,MELVA3 449 : ELSE IF (ICAS.EQ.2) THEN 450 : SEGSUP MELVA2,MELVA3 451 : SEGDES MELVA1 452 : ELSE 453 : SEGSUP MELVA1 454 : SEGDES MELVA2,MELVA3 455 : ENDIF 456 : RETURN 457 : ENDIF 458 : * 459 : * 460 : * 461 : * cas de la barre 462 : * 463 : 46 CONTINUE 464 : * ================ calcul de la vitesse du son 465 : * le resultat est stocké dans melval avec n1ptel =1 466 : IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN 467 : * recherche des paramètre matériau 468 : * module d'young 469 : MELVA3 = MPTVA1.IVAL(1) 470 : * densite 471 : MELVA4 = MPTVA1.IVAL(3) 472 : SEGDES MPTVA1 473 : SEGACT MELVA3,MELVA4 474 : * 475 : N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2)) 476 : N1PTEL = 1 477 : N2PTEL = 0 478 : N2EL = 0 479 : SEGINI MELVAL 480 : MCSON = MELVAL 481 : * boucle sur les éléments pour calculer la vitesse du son 482 : DO 4603 I=1,N1EL 483 : * on prend les valeurs moyennes sur les éléments 484 : YOU1 = 0.D0 485 : I3 = MIN(I,MELVA3.VELCHE(/2)) 486 : DO 4601 J=1,MELVA3.VELCHE(/1) 487 : YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 488 : 4601 CONTINUE 489 : YOU1 = YOU1 / MELVA3.VELCHE(/1) 490 : * 491 : RO1 = 0.D0 492 : I4 = MIN(I,MELVA4.VELCHE(/2)) 493 : DO 4602 J=1,MELVA4.VELCHE(/1) 494 : RO1 = RO1 + MELVA4.VELCHE(J,I4) 495 : 4602 CONTINUE 496 : RO1 = RO1 / MELVA4.VELCHE(/1) 497 : IF (RO1.EQ.0.D0) THEN 498 : SEGDES MELVA4,MELVA3 499 : SEGSUP MELVAL 499 : CALL ERREUR(855) 500 : RETURN 501 : ENDIF 502 : * 503 : IF (YOU1.EQ.0.D0) THEN 504 : SEGDES MELVA4,MELVA3 505 : SEGSUP MELVAL 505 : CALL ERREUR(856) 506 : RETURN 507 : ENDIF 508 : VELCHE(1,I) = SQRT(YOU1/RO1) 509 : * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 510 : 4603 CONTINUE 511 : SEGDES MELVA4,MELVA3 512 : IF (ICAS.EQ.4) THEN 513 : * cas ou seule la vitesse du son est demandée 514 : MELVA2 = 0 515 : MELV1 = MELVAL 516 : SEGDES MELVAL 517 : RETURN 518 : ENDIF 519 : ELSE IF (ICAS.EQ.2) THEN 520 : * recuperation du champ 521 : MELVA1 = MPTVA2.IVAL(1) 522 : SEGACT MELVA1 523 : MCSON = MELVA1 524 : SEGDES MPTVA2 525 : ENDIF 526 : * ================ paramètre geometrique 527 : * stocké dans un melval mtaille 528 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN 529 : MELEME = IPMAIL 530 : SEGACT MELEME 531 : N1EL = NUM(/2) 532 : N1PTEL = 1 533 : N2PTEL = 0 534 : N2EL = 0 535 : SEGINI MELVAL 536 : MTAIL1 = MELVAL 537 : MTAIL2 = 0 537 : CALL CFLTAI(MTAIL1,IPMAIL,MELE) 538 : SEGDES MELEME 539 : IF (ICAS.EQ.5) THEN 540 : MELV1 = MTAIL1 541 : MELV2 = 0 542 : SEGDES MELVAL 543 : RETURN 544 : ENDIF 545 : ELSE IF (ICAS.EQ.3) THEN 546 : * recuperation du champ 547 : MELVA1 = MPTVA2.IVAL(1) 548 : SEGACT MELVA1 549 : MTAIL1 = MELVA1 550 : SEGDES MPTVA2 551 : ENDIF 552 : * ================ pas de temps cfl 553 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN 554 : * recuperation de la vitesse du son 555 : * et du paramètre de taille 556 : MELVA1 = MCSON 557 : MELVA2 = MTAIL1 558 : * creation du melval résultat 559 : N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2)) 560 : N1PTEL = 1 561 : N2EL = 0 562 : N2PTEL = 0 563 : SEGINI MELVAL 564 : * 565 : DO 4604 I=1,N1EL 566 : I1 = MIN(I,MELVA1.VELCHE(/2)) 567 : I2 = MIN(I,MELVA2.VELCHE(/2)) 568 : VELCHE(1,I)=MELVA2.VELCHE(1,I2)/MELVA1.VELCHE(1,I1) 569 : * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 570 : 4604 CONTINUE 571 : MELV1 = MELVAL 572 : MELV2 = 0 573 : SEGDES MELVAL 574 : IF (ICAS.EQ.1) THEN 575 : SEGSUP MELVA1,MELVA2 576 : ELSE IF (ICAS.EQ.2) THEN 577 : SEGSUP MELVA2 578 : SEGDES MELVA1 579 : ELSE 580 : SEGSUP MELVA1 581 : SEGDES MELVA2 582 : ENDIF 583 : RETURN 584 : ENDIF 585 : * 586 : * 587 : * cas du coq3 et du dkt (meme avec excentrement) 588 : * Rq: dst aussi dans les cas simples 589 : * 590 : 27 CONTINUE 591 : * ================ calcul de la vitesse du son 592 : * le resultat est stocké dans melval avec n1ptel =1 593 : IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN 594 : * recherche des paramètre matériau 595 : * module d'young 596 : MELVA3 = MPTVA1.IVAL(1) 597 : * densite 598 : MELVA4 = MPTVA1.IVAL(3) 599 : * nu 600 : MELVA5 = MPTVA1.IVAL(2) 601 : * SEGDES MPTVA1 602 : SEGACT MELVA3,MELVA4,MELVA5 603 : * 604 : N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2), 605 : & MELVA5.VELCHE(/2)) 606 : N1PTEL = 1 607 : N2PTEL = 0 608 : N2EL = 0 609 : SEGINI MELVAL 610 : MCSON = MELVAL 611 : * boucle sur les éléments pour calculer la vitesse du son 612 : DO 2704 I=1,N1EL 613 : * on prend les valeurs moyennes sur les éléments 614 : YOU1 = 0.D0 615 : I3 = MIN(I,MELVA3.VELCHE(/2)) 616 : DO 2701 J=1,MELVA3.VELCHE(/1) 617 : YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 618 : 2701 CONTINUE 619 : YOU1 = YOU1 / MELVA3.VELCHE(/1) 620 : * 621 : RO1 = 0.D0 622 : I4 = MIN(I,MELVA4.VELCHE(/2)) 623 : DO 2702 J=1,MELVA4.VELCHE(/1) 624 : RO1 = RO1 + MELVA4.VELCHE(J,I4) 625 : 2702 CONTINUE 626 : RO1 = RO1 / MELVA4.VELCHE(/1) 627 : * 628 : DNU1 = 0.D0 629 : I5 = MIN(I,MELVA5.VELCHE(/2)) 630 : DO 2703 J=1,MELVA5.VELCHE(/1) 631 : DNU1 = DNU1 + MELVA5.VELCHE(J,I5) 632 : 2703 CONTINUE 633 : DNU1 = DNU1 / MELVA5.VELCHE(/1) 634 : 635 : IF (RO1.EQ.0.D0) THEN 636 : SEGDES MELVA4,MELVA3,MELVA5 637 : SEGSUP MELVAL 637 : CALL ERREUR(855) 638 : RETURN 639 : ENDIF 640 : * 641 : IF (YOU1.EQ.0.D0) THEN 642 : SEGDES MELVA4,MELVA3,MELVA5 643 : SEGSUP MELVAL 643 : CALL ERREUR(856) 644 : RETURN 645 : ENDIF 646 : IF (DNU1.GE.1.D0) THEN 647 : SEGDES MELVA4,MELVA3,MELVA5 648 : SEGSUP MELVAL 648 : CALL ERREUR(36) 649 : RETURN 650 : ENDIF 651 : * 652 : VELCHE(1,I) = SQRT(YOU1/(RO1*(1-DNU1*DNU1))) 653 : * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 654 : 2704 CONTINUE 655 : SEGDES MELVA4,MELVA3,MELVA5 656 : IF (ICAS.EQ.4) THEN 657 : * cas ou seule la vitesse du son est demandée 658 : MELVA2 = 0 659 : MELV1 = MELVAL 660 : SEGDES MELVAL 661 : RETURN 662 : ENDIF 663 : ELSE IF (ICAS.EQ.2) THEN 664 : * recuperation du champ 665 : MELVA1 = MPTVA2.IVAL(1) 666 : SEGACT MELVA1 667 : MCSON = MELVA1 668 : SEGDES MPTVA2 669 : ENDIF 670 : * ================ paramètre geometrique 671 : * stocké dans un melval mtaille 672 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN 673 : MELEME = IPMAIL 674 : SEGACT MELEME 675 : N1EL = NUM(/2) 676 : N1PTEL = 1 677 : N2PTEL = 0 678 : N2EL = 0 679 : SEGINI MELVAL 680 : MTAIL1 = MELVAL 680 : CALL CFLTAI(MTAIL1,IPMAIL,MELE) 681 : SEGDES MELEME 682 : * second paramètre geométrique 683 : IOBL = 3 684 : IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0 685 : SEGINI, MELVA2=MELVAL 686 : MTAIL2 = MELVA2 687 : * récupération de l' epaisseur 688 : MELVA3 = MPTVA1.IVAL(1+IOBL) 689 : SEGACT MELVA3 690 : SEGACT MELEME 691 : DO 2706 I=1,N1EL 692 : I3 = MIN(I,MELVA3.VELCHE(/2)) 693 : * remarque les valeurs ne sont pas moyennées sur l'élément 694 : V12L = 0.D0 695 : V13L = 0.D0 696 : V23L = 0.D0 697 : DO 2761 J=1,IDIM 698 : V13(J)=XCOOR((IDIM+1)*(NUM(3,I3)-1)+J)- 699 : & XCOOR((IDIM+1)*(NUM(1,I3)-1)+J) 700 : V12(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)- 701 : & XCOOR((IDIM+1)*(NUM(1,I)-1)+J) 702 : V23(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)- 703 : & XCOOR((IDIM+1)*(NUM(2,I)-1)+J) 704 : V12L = V12L + V12(J)*V12(J) 705 : V13L = V13L + V13(J)*V13(J) 706 : V23L = V23L + V23(J)*V23(J) 707 : 2761 CONTINUE 708 : D2min = 0.D0 709 : D2max = 0.D0 710 : Dmax = 0.D0 711 : D2min = MIN(V12L,V13L,V23L) 712 : D2max = MAX(V12L,V13L,V23L) 713 : Dmax = SQRT(D2max) 714 : * write(6,*) 'D2min', D2min , 'D2max' , D2max, 'Dmax', Dmax 715 : * write(6,*) 'Element', i , 'L' , MELVA2.VELCHE(1,i) 716 : DUM1 = SQRT(D2min-MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I)) 717 : DUM1 = DUM1 / Dmax - 0.5D0 718 : DUM1 = 5.3706D0 * DUM1 * DUM1 719 : DUM1 = 3.9402D0 - 0.8687D0*MELVA2.VELCHE(1,I)/Dmax - DUM1 720 : DUM1 = Dmax / (XPI * MELVA3.VELCHE(1,I3) * DUM1) 721 : DUM1 = DUM1 * MELVA2.VELCHE(1,I) 722 : IF (MELE.EQ.28 .AND. MPTVA1.IVAL(2+IOBL).NE.0) THEN 723 : * récupération de l'excentricité 724 : MELVA4 = MPTVA1.IVAL(2+IOBL) 725 : SEGACT MELVA4 726 : I4 = MIN(I,MELVA4.VELCHE(/2)) 727 : DUM1=DUM1/SQRT(1D0+12D0*MELVA4.VELCHE(1,I4)*MELVA4.VELCHE(1,I4) 728 : & /(MELVA3.VELCHE(1,I3)*MELVA3.VELCHE(1,I3))) 729 : * On applique un coefficient empirique de 1.5 qui tient compte des 730 : * différences de discrétisation par rapport au modèle non excentré 731 : DUM1 = DUM1 * 1.5D0 732 : SEGDES MELVA4 733 : ENDIF 734 : MELVA2.VELCHE(1,I)=DUM1 735 : 2706 CONTINUE 736 : SEGDES MELEME 737 : SEGDES MELVA3 738 : SEGDES MPTVA1 739 : * 740 : IF (ICAS.EQ.5) THEN 741 : MELV1 = MTAIL1 742 : MELV2 = MTAIL2 743 : SEGDES MELVAL 744 : RETURN 745 : ENDIF 746 : ELSE IF (ICAS.EQ.3) THEN 747 : * recuperation du champ 748 : MELVA1 = MPTVA2.IVAL(1) 749 : MELVA2 = MPTVA2.IVAL(2) 750 : SEGACT MELVA1,MELVA2 751 : MTAIL1 = MELVA1 752 : MTAIL2 = MELVA2 753 : SEGDES MPTVA2 754 : ENDIF 755 : * ================ pas de temps cfl 756 : IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN 757 : * recuperation de la vitesse du son 758 : * et du paramètre de taille 759 : MELVA1 = MCSON 760 : MELVA2 = MTAIL1 761 : MELVA3 = MTAIL2 762 : * creation du melval résultat 763 : N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2)) 764 : N1PTEL = 1 765 : N2EL = 0 766 : N2PTEL = 0 767 : SEGINI MELVAL 768 : * 769 : DO 2705 I=1,N1EL 770 : I1 = MIN(I,MELVA1.VELCHE(/2)) 771 : I2 = MIN(I,MELVA2.VELCHE(/2)) 772 : I3 = MIN(I,MELVA3.VELCHE(/2)) 773 : DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3)) 774 : VELCHE(1,I)=DUM1/MELVA1.VELCHE(1,I1) 775 : * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 776 : 2705 CONTINUE 777 : MELV1 = MELVAL 778 : MELV2 = 0 779 : SEGDES MELVAL 780 : IF (ICAS.EQ.1) THEN 781 : SEGSUP MELVA1,MELVA2,MELVA3 782 : ELSE IF (ICAS.EQ.2) THEN 783 : SEGSUP MELVA2,MELVA3 784 : SEGDES MELVA1 785 : ELSE 786 : SEGSUP MELVA1 787 : SEGDES MELVA2,MELVA3 788 : ENDIF 789 : RETURN 790 : ENDIF 791 : * 792 : 99 CONTINUE 793 : MOTERR(1:4)=NOMTP(MELE) 794 : MOTERR(9:12)='CFL5' 794 : CALL ERREUR(86) 795 : * 796 : RETURN 797 : END 798 : 799 : 800 : 801 : 802 : 803 : 804 : 805 : 806 :
© Cast3M 2003 - All rights reserved.
Disclaimer