Télécharger la plaquette Cast3M


fr en

Télécharger cfl5.eso

Retour à la liste

Numérotation des lignes : 

cfl5

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales