C CFL5 SOURCE OF166741 25/02/21 21:15:30 12166 SUBROUTINE CFL5(ICAS,IPMAIL,MELE,IVAM1,IVAM2,MELV1,MELV2,N2) *---------------------------------------------------------------------* * * calcul du pas de temps CFL * * elements poutre,tuyau,barre,cerc,coque2,coque3,dkt * * * entree * icas : cas à traiter * = 1 calcul du pas de temps complet ivam1 avec cara * = 2 calcul du pas de temps lorsque cson est donne ivam2 * = 3 calcul du pas de temps lorsque la taille est donnée ivam1 si cara * = 4 calcul de la vitesse du son ivam1 donné * = 5 calcul du parametre de taille ivam1 si cara * * ipmail : pointeur vers le maillage a traiter * mele : numero de l'élément finis dans nomtp * ivam1 : pointeur vers mptval du cham1 actif * ivam2 : pointeur vers mptval du cham2 actif * n2 : nombre de composante en sortie * * sortie * melv1 : melval de la première composante du chamelem resultat * inactif en sortie * melv2 : melval de la deuxième composante du chamelem resultat * inactif en sortie * * *---------------------------------------------------------------------* IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC CCHAMP -INC SMCHAML -INC SMELEME -INC SMCOORD -INC TMPTVAL DIMENSION V13(3),V12(3),V23(3) MPTVA1 = IVAM1 MPTVA2 = IVAM2 * * * branchement en fonction de l'élément fini * * 0 5 0 5 0 GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 2 99,99,99,99,99,99,27,27,29,99,99,99,99,99,99,99,99,99,99,99, * bar 4 99,29,99,44,99,46,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 6 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 8 99,99,99,29,99,99,99,99,99,99,99,99,27,99,99,99,99,99,99,99, 1 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 2 99,99,99,99,99,99,99),MELE GOTO 99 * * cas de la poutre POUT (mince) et de TIMO (épaisse) * 29 CONTINUE * ================ calcul de la vitesse du son * le resultat est stocké dans melval avec n1ptel =1 IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN * recherche des paramètre matériau * module d'young MELVA3 = MPTVA1.IVAL(1) * densite MELVA4 = MPTVA1.IVAL(3) SEGACT MELVA3,MELVA4 * N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2)) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MCSON = MELVAL * boucle sur les éléments pour calculer la vitesse du son DO 2903 I=1,N1EL * on prend les valeurs moyennes sur les éléments YOU1 = 0.D0 I3 = MIN(I,MELVA3.VELCHE(/2)) DO 2901 J=1,MELVA3.VELCHE(/1) YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 2901 CONTINUE YOU1 = YOU1 / MELVA3.VELCHE(/1) * RO1 = 0.D0 I4 = MIN(I,MELVA4.VELCHE(/2)) DO 2902 J=1,MELVA4.VELCHE(/1) RO1 = RO1 + MELVA4.VELCHE(J,I4) 2902 CONTINUE RO1 = RO1 / MELVA4.VELCHE(/1) IF (RO1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3 SEGSUP MELVAL CALL ERREUR(855) RETURN ENDIF * IF (YOU1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3 SEGSUP MELVAL CALL ERREUR(856) RETURN ENDIF VELCHE(1,I) = SQRT(YOU1/RO1) * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 2903 CONTINUE SEGDES MELVA4,MELVA3 IF (ICAS.EQ.4) THEN * cas ou seule la vitesse du son est demandée MELV2 = 0 MELV1 = MELVAL SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.2) THEN * recuperation du champ * SEGACT MPTVA2 MELVA1 = MPTVA2.IVAL(1) SEGACT MELVA1 MCSON = MELVA1 SEGDES MPTVA2 ENDIF * ================ paramètre geometrique * stocké dans un melval mtaille IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN * on verifie que le champ de caractéristiques existe IF (IVAM1.EQ.0) THEN MOTERR(1:8)='CARACTER' CALL ERREUR(145) RETURN ENDIF MELEME = IPMAIL SEGACT MELEME N1EL = NUM(/2) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MTAIL1 = MELVAL CALL CFLTAI(MTAIL1,IPMAIL,MELE) SEGDES MELEME * second paramètre geométrique IOBL = 3 IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0 SEGINI, MELVA2=MELVAL MTAIL2 = MELVA2 IF (MELE.EQ.29) THEN * récupération de la section des inerties et de la torsion * section MELVA3 = MPTVA1.IVAL(4+IOBL) * inertie z MELVA4 = MPTVA1.IVAL(3+IOBL) * inertie y MELVA5 = MPTVA1.IVAL(2+IOBL) * torsion MELVA6 = MPTVA1.IVAL(1+IOBL) SEGACT MELVA3,MELVA4,MELVA5,MELVA6 DO 2905 I=1,N1EL I3 = MIN(I,MELVA3.VELCHE(/2)) I4 = MIN(I,MELVA4.VELCHE(/2)) I5 = MIN(I,MELVA5.VELCHE(/2)) I6 = MIN(I,MELVA6.VELCHE(/2)) * remaarque les valeurs ne sont pas moyennées sur l'élément DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) DUM1 = DUM1 * MELVA3.VELCHE(1,I3) DUM1=DUM1/MAX(MELVA4.VELCHE(1,I4),MELVA5.VELCHE(1,I5), & MELVA6.VELCHE(1,I6)) DUM1=SQRT(DUM1) * 0.10132D0 MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 2905 CONTINUE SEGDES MELVA3,MELVA4,MELVA5,MELVA6 ELSEIF (MELE.EQ.84) THEN * récupération de la section des inerties et de la torsion * section MELVA3 = MPTVA1.IVAL(4+IOBL) * inertie z MELVA4 = MPTVA1.IVAL(3+IOBL) * inertie y MELVA5 = MPTVA1.IVAL(2+IOBL) * torsion MELVA6 = MPTVA1.IVAL(1+IOBL) SEGACT MELVA3,MELVA4,MELVA5,MELVA6 DO 8405 I=1,N1EL I3 = MIN(I,MELVA3.VELCHE(/2)) I4 = MIN(I,MELVA4.VELCHE(/2)) I5 = MIN(I,MELVA5.VELCHE(/2)) I6 = MIN(I,MELVA6.VELCHE(/2)) * remarque les valeurs ne sont pas moyennées sur l'élément DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) * 0.9 = 2**1.5 / XPI 0.7 coef de securite DUM2 = MELVA2.VELCHE(1,I)*0.9D0*0.7D0 DUM1 = DUM1 * MELVA3.VELCHE(1,I3) DUM1=DUM1/MAX(MELVA4.VELCHE(1,I4),MELVA5.VELCHE(1,I5), & MELVA6.VELCHE(1,I6)) DUM1=SQRT(DUM1) * 0.10132D0 * correction pour tenir compte du cisaillement Blevins page 175 eq 8.30 MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 + DUM2 8405 CONTINUE SEGDES MELVA3,MELVA4,MELVA5,MELVA6 ELSEIF(MELE.EQ.42) THEN * cas du tuyau * récupération des caractéristiques * epai MELVA3 = MPTVA1.IVAL(1+IOBL) * rayo MELVA4 = MPTVA1.IVAL(2+IOBL) * pour les autres on ne ient pas compte des modification * qui assouplissent le tuyau donc omega max diminue SEGACT MELVA3,MELVA4 DO 4205 I=1,N1EL I3 = MIN(I,MELVA3.VELCHE(/2)) I4 = MIN(I,MELVA4.VELCHE(/2)) REXT = MELVA4.VELCHE(1,I4) RINT = REXT - MELVA3.VELCHE(1,I3) * remarque les valeurs ne sont pas moyennées sur l'élément DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) DUM1 = DUM1 * 2.D0 / (REXT*REXT+RINT*RINT) DUM1=SQRT(DUM1) * 0.10132D0 MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 4205 CONTINUE SEGDES MELVA3,MELVA4 ENDIF IF (ICAS.EQ.5) THEN MELV1 = MTAIL1 MELV2 = MTAIL2 SEGDES MELVAL,MELVA2 RETURN ENDIF ELSE IF (ICAS.EQ.3) THEN * recuperation du champ * SEGACT MPTVA2 MELVA1 = MPTVA2.IVAL(1) MELVA2 = MPTVA2.IVAL(2) SEGACT MELVA1,MELVA2 MTAIL1 = MELVA1 MTAIL2 = MELVA2 SEGDES MPTVA2 ENDIF * ================ pas de temps cfl IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN * recuperation de la vitesse du son * et du paramètre de taille MELVA1 = MCSON MELVA2 = MTAIL1 MELVA3 = MTAIL2 * creation du melval résultat N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2), & MELVA3.VELCHE(/2)) N1PTEL = 1 N2EL = 0 N2PTEL = 0 SEGINI MELVAL * DO 2904 I=1,N1EL I1 = MIN(I,MELVA1.VELCHE(/2)) I2 = MIN(I,MELVA2.VELCHE(/2)) I3 = MIN(I,MELVA3.VELCHE(/2)) DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3)) VELCHE(1,I)= DUM1 / MELVA1.VELCHE(1,I1) * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 2904 CONTINUE MELV1 = MELVAL MELV2 = 0 SEGDES MELVAL IF (ICAS.EQ.1) THEN SEGSUP MELVA1,MELVA2,MELVA3 ELSE IF (ICAS.EQ.2) THEN SEGSUP MELVA2,MELVA3 SEGDES MELVA1 ELSE * icas=3 SEGSUP MELVA1 SEGDES MELVA2,MELVA3 ENDIF RETURN ENDIF * * * cas du coq2 * 44 CONTINUE * ================ calcul de la vitesse du son * le resultat est stocké dans melval avec n1ptel =1 IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN * recherche des paramètre matériau * module d'young MELVA3 = MPTVA1.IVAL(1) * densite MELVA4 = MPTVA1.IVAL(3) * nu MELVA5 = MPTVA1.IVAL(2) * SEGDES MPTVA1 SEGACT MELVA3,MELVA4,MELVA5 * N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2), & MELVA5.VELCHE(/2)) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MCSON = MELVAL * boucle sur les éléments pour calculer la vitesse du son DO 4404 I=1,N1EL * on prend les valeurs moyennes sur les éléments YOU1 = 0.D0 I3 = MIN(I,MELVA3.VELCHE(/2)) DO 4401 J=1,MELVA3.VELCHE(/1) YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 4401 CONTINUE YOU1 = YOU1 / MELVA3.VELCHE(/1) * RO1 = 0.D0 I4 = MIN(I,MELVA4.VELCHE(/2)) DO 4402 J=1,MELVA4.VELCHE(/1) RO1 = RO1 + MELVA4.VELCHE(J,I4) 4402 CONTINUE RO1 = RO1 / MELVA4.VELCHE(/1) * DNU1 = 0.D0 I5 = MIN(I,MELVA5.VELCHE(/2)) DO 4403 J=1,MELVA5.VELCHE(/1) DNU1 = DNU1 + MELVA5.VELCHE(J,I5) 4403 CONTINUE DNU1 = DNU1 / MELVA5.VELCHE(/1) IF (RO1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3,MELVA5 SEGSUP MELVAL CALL ERREUR(855) RETURN ENDIF * IF (YOU1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3,MELVA5 SEGSUP MELVAL CALL ERREUR(856) RETURN ENDIF IF (DNU1.GE.1.D0) THEN SEGDES MELVA4,MELVA3,MELVA5 SEGSUP MELVAL CALL ERREUR(36) RETURN ENDIF * VELCHE(1,I) = SQRT(YOU1/(RO1*(1-DNU1*DNU1))) * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 4404 CONTINUE SEGDES MELVA4,MELVA3,MELVA5 IF (ICAS.EQ.4) THEN * cas ou seule la vitesse du son est demandée MELVA2 = 0 MELV1 = MELVAL SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.2) THEN * recuperation du champ MELVA1 = MPTVA2.IVAL(1) SEGACT MELVA1 MCSON = MELVA1 SEGDES MPTVA2 ENDIF * ================ paramètre geometrique * stocké dans un melval mtaille IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN MELEME = IPMAIL SEGACT MELEME N1EL = NUM(/2) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MTAIL1 = MELVAL CALL CFLTAI(MTAIL1,IPMAIL,MELE) SEGDES MELEME * second paramètre geométrique IOBL = 3 IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0 SEGINI, MELVA2=MELVAL MTAIL2 = MELVA2 * récupération de l' epaisseur MELVA3 = MPTVA1.IVAL(1+IOBL) SEGACT MELVA3 DO 4406 I=1,N1EL I3 = MIN(I,MELVA3.VELCHE(/2)) * remarque les valeurs ne sont pas moyennées sur l'élément DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I) DUM1 = DUM1 / (3.23D0 * MELVA3.VELCHE(1,I3)) MELVA2.VELCHE(1,I)=DUM1 4406 CONTINUE SEGDES MELVA3 SEGDES MPTVA1 * IF (ICAS.EQ.5) THEN MELV1 = MTAIL1 MELV2 = MTAIL2 SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.3) THEN * recuperation du champ MELVA1 = MPTVA2.IVAL(1) MELVA2 = MPTVA2.IVAL(2) SEGACT MELVA1,MELVA2 MTAIL1 = MELVA1 MTAIL2 = MELVA2 SEGDES MPTVA2 ENDIF * ================ pas de temps cfl IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN * recuperation de la vitesse du son * et du paramètre de taille MELVA1 = MCSON MELVA2 = MTAIL1 MELVA3 = MTAIL2 * creation du melval résultat N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2)) N1PTEL = 1 N2EL = 0 N2PTEL = 0 SEGINI MELVAL * DO 4405 I=1,N1EL I1 = MIN(I,MELVA1.VELCHE(/2)) I2 = MIN(I,MELVA2.VELCHE(/2)) I3 = MIN(I,MELVA3.VELCHE(/2)) DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3)) VELCHE(1,I)=DUM1/MELVA1.VELCHE(1,I1) * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 4405 CONTINUE MELV1 = MELVAL MELV2 = 0 SEGDES MELVAL IF (ICAS.EQ.1) THEN SEGSUP MELVA1,MELVA2,MELVA3 ELSE IF (ICAS.EQ.2) THEN SEGSUP MELVA2,MELVA3 SEGDES MELVA1 ELSE SEGSUP MELVA1 SEGDES MELVA2,MELVA3 ENDIF RETURN ENDIF * * * * cas de la barre * 46 CONTINUE * ================ calcul de la vitesse du son * le resultat est stocké dans melval avec n1ptel =1 IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN * recherche des paramètre matériau * module d'young MELVA3 = MPTVA1.IVAL(1) * densite MELVA4 = MPTVA1.IVAL(3) SEGDES MPTVA1 SEGACT MELVA3,MELVA4 * N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2)) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MCSON = MELVAL * boucle sur les éléments pour calculer la vitesse du son DO 4603 I=1,N1EL * on prend les valeurs moyennes sur les éléments YOU1 = 0.D0 I3 = MIN(I,MELVA3.VELCHE(/2)) DO 4601 J=1,MELVA3.VELCHE(/1) YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 4601 CONTINUE YOU1 = YOU1 / MELVA3.VELCHE(/1) * RO1 = 0.D0 I4 = MIN(I,MELVA4.VELCHE(/2)) DO 4602 J=1,MELVA4.VELCHE(/1) RO1 = RO1 + MELVA4.VELCHE(J,I4) 4602 CONTINUE RO1 = RO1 / MELVA4.VELCHE(/1) IF (RO1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3 SEGSUP MELVAL CALL ERREUR(855) RETURN ENDIF * IF (YOU1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3 SEGSUP MELVAL CALL ERREUR(856) RETURN ENDIF VELCHE(1,I) = SQRT(YOU1/RO1) * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 4603 CONTINUE SEGDES MELVA4,MELVA3 IF (ICAS.EQ.4) THEN * cas ou seule la vitesse du son est demandée MELVA2 = 0 MELV1 = MELVAL SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.2) THEN * recuperation du champ MELVA1 = MPTVA2.IVAL(1) SEGACT MELVA1 MCSON = MELVA1 SEGDES MPTVA2 ENDIF * ================ paramètre geometrique * stocké dans un melval mtaille IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN MELEME = IPMAIL SEGACT MELEME N1EL = NUM(/2) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MTAIL1 = MELVAL MTAIL2 = 0 CALL CFLTAI(MTAIL1,IPMAIL,MELE) SEGDES MELEME IF (ICAS.EQ.5) THEN MELV1 = MTAIL1 MELV2 = 0 SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.3) THEN * recuperation du champ MELVA1 = MPTVA2.IVAL(1) SEGACT MELVA1 MTAIL1 = MELVA1 SEGDES MPTVA2 ENDIF * ================ pas de temps cfl IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN * recuperation de la vitesse du son * et du paramètre de taille MELVA1 = MCSON MELVA2 = MTAIL1 * creation du melval résultat N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2)) N1PTEL = 1 N2EL = 0 N2PTEL = 0 SEGINI MELVAL * DO 4604 I=1,N1EL I1 = MIN(I,MELVA1.VELCHE(/2)) I2 = MIN(I,MELVA2.VELCHE(/2)) VELCHE(1,I)=MELVA2.VELCHE(1,I2)/MELVA1.VELCHE(1,I1) * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 4604 CONTINUE MELV1 = MELVAL MELV2 = 0 SEGDES MELVAL IF (ICAS.EQ.1) THEN SEGSUP MELVA1,MELVA2 ELSE IF (ICAS.EQ.2) THEN SEGSUP MELVA2 SEGDES MELVA1 ELSE SEGSUP MELVA1 SEGDES MELVA2 ENDIF RETURN ENDIF * * * cas du coq3 et du dkt (meme avec excentrement) * Rq: dst aussi dans les cas simples * 27 CONTINUE * ================ calcul de la vitesse du son * le resultat est stocké dans melval avec n1ptel =1 IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN * recherche des paramètre matériau * module d'young MELVA3 = MPTVA1.IVAL(1) * densite MELVA4 = MPTVA1.IVAL(3) * nu MELVA5 = MPTVA1.IVAL(2) * SEGDES MPTVA1 SEGACT MELVA3,MELVA4,MELVA5 * N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2), & MELVA5.VELCHE(/2)) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MCSON = MELVAL * boucle sur les éléments pour calculer la vitesse du son DO 2704 I=1,N1EL * on prend les valeurs moyennes sur les éléments YOU1 = 0.D0 I3 = MIN(I,MELVA3.VELCHE(/2)) DO 2701 J=1,MELVA3.VELCHE(/1) YOU1 = YOU1 + MELVA3.VELCHE(J,I3) 2701 CONTINUE YOU1 = YOU1 / MELVA3.VELCHE(/1) * RO1 = 0.D0 I4 = MIN(I,MELVA4.VELCHE(/2)) DO 2702 J=1,MELVA4.VELCHE(/1) RO1 = RO1 + MELVA4.VELCHE(J,I4) 2702 CONTINUE RO1 = RO1 / MELVA4.VELCHE(/1) * DNU1 = 0.D0 I5 = MIN(I,MELVA5.VELCHE(/2)) DO 2703 J=1,MELVA5.VELCHE(/1) DNU1 = DNU1 + MELVA5.VELCHE(J,I5) 2703 CONTINUE DNU1 = DNU1 / MELVA5.VELCHE(/1) IF (RO1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3,MELVA5 SEGSUP MELVAL CALL ERREUR(855) RETURN ENDIF * IF (YOU1.EQ.0.D0) THEN SEGDES MELVA4,MELVA3,MELVA5 SEGSUP MELVAL CALL ERREUR(856) RETURN ENDIF IF (DNU1.GE.1.D0) THEN SEGDES MELVA4,MELVA3,MELVA5 SEGSUP MELVAL CALL ERREUR(36) RETURN ENDIF * VELCHE(1,I) = SQRT(YOU1/(RO1*(1-DNU1*DNU1))) * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i) 2704 CONTINUE SEGDES MELVA4,MELVA3,MELVA5 IF (ICAS.EQ.4) THEN * cas ou seule la vitesse du son est demandée MELVA2 = 0 MELV1 = MELVAL SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.2) THEN * recuperation du champ MELVA1 = MPTVA2.IVAL(1) SEGACT MELVA1 MCSON = MELVA1 SEGDES MPTVA2 ENDIF * ================ paramètre geometrique * stocké dans un melval mtaille IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN MELEME = IPMAIL SEGACT MELEME N1EL = NUM(/2) N1PTEL = 1 N2PTEL = 0 N2EL = 0 SEGINI MELVAL MTAIL1 = MELVAL CALL CFLTAI(MTAIL1,IPMAIL,MELE) SEGDES MELEME * second paramètre geométrique IOBL = 3 IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0 SEGINI, MELVA2=MELVAL MTAIL2 = MELVA2 * récupération de l' epaisseur MELVA3 = MPTVA1.IVAL(1+IOBL) SEGACT MELVA3 SEGACT MELEME DO 2706 I=1,N1EL I3 = MIN(I,MELVA3.VELCHE(/2)) * remarque les valeurs ne sont pas moyennées sur l'élément V12L = 0.D0 V13L = 0.D0 V23L = 0.D0 DO 2761 J=1,IDIM V13(J)=XCOOR((IDIM+1)*(NUM(3,I3)-1)+J)- & XCOOR((IDIM+1)*(NUM(1,I3)-1)+J) V12(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)- & XCOOR((IDIM+1)*(NUM(1,I)-1)+J) V23(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)- & XCOOR((IDIM+1)*(NUM(2,I)-1)+J) V12L = V12L + V12(J)*V12(J) V13L = V13L + V13(J)*V13(J) V23L = V23L + V23(J)*V23(J) 2761 CONTINUE D2min = 0.D0 D2max = 0.D0 Dmax = 0.D0 D2min = MIN(V12L,V13L,V23L) D2max = MAX(V12L,V13L,V23L) Dmax = SQRT(D2max) * write(6,*) 'D2min', D2min , 'D2max' , D2max, 'Dmax', Dmax * write(6,*) 'Element', i , 'L' , MELVA2.VELCHE(1,i) DUM1 = SQRT(D2min-MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I)) DUM1 = DUM1 / Dmax - 0.5D0 DUM1 = 5.3706D0 * DUM1 * DUM1 DUM1 = 3.9402D0 - 0.8687D0*MELVA2.VELCHE(1,I)/Dmax - DUM1 DUM1 = Dmax / (XPI * MELVA3.VELCHE(1,I3) * DUM1) DUM1 = DUM1 * MELVA2.VELCHE(1,I) IF (MELE.EQ.28 .AND. MPTVA1.IVAL(2+IOBL).NE.0) THEN * récupération de l'excentricité MELVA4 = MPTVA1.IVAL(2+IOBL) SEGACT MELVA4 I4 = MIN(I,MELVA4.VELCHE(/2)) DUM1=DUM1/SQRT(1D0+12D0*MELVA4.VELCHE(1,I4)*MELVA4.VELCHE(1,I4) & /(MELVA3.VELCHE(1,I3)*MELVA3.VELCHE(1,I3))) * On applique un coefficient empirique de 1.5 qui tient compte des * différences de discrétisation par rapport au modèle non excentré DUM1 = DUM1 * 1.5D0 SEGDES MELVA4 ENDIF MELVA2.VELCHE(1,I)=DUM1 2706 CONTINUE SEGDES MELEME SEGDES MELVA3 SEGDES MPTVA1 * IF (ICAS.EQ.5) THEN MELV1 = MTAIL1 MELV2 = MTAIL2 SEGDES MELVAL RETURN ENDIF ELSE IF (ICAS.EQ.3) THEN * recuperation du champ MELVA1 = MPTVA2.IVAL(1) MELVA2 = MPTVA2.IVAL(2) SEGACT MELVA1,MELVA2 MTAIL1 = MELVA1 MTAIL2 = MELVA2 SEGDES MPTVA2 ENDIF * ================ pas de temps cfl IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN * recuperation de la vitesse du son * et du paramètre de taille MELVA1 = MCSON MELVA2 = MTAIL1 MELVA3 = MTAIL2 * creation du melval résultat N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2)) N1PTEL = 1 N2EL = 0 N2PTEL = 0 SEGINI MELVAL * DO 2705 I=1,N1EL I1 = MIN(I,MELVA1.VELCHE(/2)) I2 = MIN(I,MELVA2.VELCHE(/2)) I3 = MIN(I,MELVA3.VELCHE(/2)) DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3)) VELCHE(1,I)=DUM1/MELVA1.VELCHE(1,I1) * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i) 2705 CONTINUE MELV1 = MELVAL MELV2 = 0 SEGDES MELVAL IF (ICAS.EQ.1) THEN SEGSUP MELVA1,MELVA2,MELVA3 ELSE IF (ICAS.EQ.2) THEN SEGSUP MELVA2,MELVA3 SEGDES MELVA1 ELSE SEGSUP MELVA1 SEGDES MELVA2,MELVA3 ENDIF RETURN ENDIF * 99 CONTINUE MOTERR(1:4)=NOMTP(MELE) MOTERR(9:12)='CFL5' CALL ERREUR(86) * RETURN END