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

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