Télécharger la plaquette Cast3M


fr en

Télécharger devfb3.eso

Retour à la liste

Numérotation des lignes : 

devfb3

  1. C DEVFB3    SOURCE    BP208322  20/09/18    21:15:28     10718          
  2.       SUBROUTINE DEVFB3(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB,
  3.      &                  NPLB,IND,IND1,INDM1,NPA,NPAM1,IND2,PDT,PDTS2,
  4.      &                  FEXPSM,NPC1,I,iannul,
  5.      &                  KTOTXB,KTOTVB,IDIMB,GETJAC)
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8(A-H,O-Z)
  8. *--------------------------------------------------------------------*
  9. *                                                                    *
  10. *     Operateur DYNE et DYNC                                         *
  11. *     Calcul des forces de choc pour les liaisons B de type :        *
  12. *     - POINT_CERCLE avec ou sans amortissement                      *
  13. *     - POINT_CERCLE_FROTTEMENT avec ou sans amortissement           *
  14. *     - POINT_CERCLE_MOBILE     avec ou sans amortissement           *
  15. *     - CERCLE_CERCLE_FROTTEMENT avec ou sans amortissement          *
  16. *                                                                    *
  17. *--------------------------------------------------------------------*
  18. *                                                                    *
  19. *     Parametres:                                                    *
  20. *                                                                    *
  21. * e   ityp    type de la liaison.                                    *
  22. * es  ftotb   forces exterieures totalisees sur la base b.           *
  23. * e   xptb    tableau des deplacements des points                    *
  24. * e   ipalb   renseigne sur la liaison.                              *
  25. * e   iplib   tableau contenant les numeros "dyne" de la liaison.    *
  26. * e   xpalb   tableau contenant les parametres de la liaison.        *
  27. * es  xvalb   tableau contenant les variables internes de liaisons.  *
  28. * e   nliab   nombre de liaisons sur la base b.                      *
  29. * e   nplb    nombre total de points intervenant dans les liaisons.  *
  30. * e   ind     indice du pas.                                         *
  31. * e   i       numero de la liaison.                                  *
  32. *                                                                    *
  33. *--------------------------------------------------------------------*
  34. *
  35. -INC CCREEL
  36.       INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*)
  37.       REAL*8  XPALB(NLIAB,*),XPTB(NPLB,2,*),FTOTB(NPLB,*)
  38.       REAL*8  XVALB(NLIAB,4,*),FEXPSM(NPLB,NPC1,2,*)
  39.       REAL*8  XPTP2(3),XPTPM2(3),XFNT(3)
  40.       REAL*8  XVPC0(3),XN2(3),XVPCT(3)
  41. *     en + pour DYNC :      
  42.       REAL*8  KTOTXB(NPLB,IDIMB,IDIMB),KTOTVB(NPLB,IDIMB,IDIMB)  
  43.       LOGICAL GETJAC
  44.  
  45.  
  46. *--------------------------------------------------------------------*
  47. * --- choc elementaire point_cercle
  48. *                      avec ou sans amortissement
  49. *--------------------------------------------------------------------*
  50. *
  51.       IF (ITYP.EQ.21 .OR. ITYP.EQ.22) THEN
  52.  
  53.          NPOI  = IPLIB(I,1)
  54.          IDIM  = IPALB(I,3)
  55.          XRAID = XPALB(I,1)
  56.          XJEU  = XPALB(I,2)
  57.  
  58.          IF (ITYP.EQ.21) THEN
  59.             ID1 = 2
  60.          ELSE
  61.             ID1 = 3
  62.          ENDIF
  63.          ID2 = ID1 + IDIM
  64. c        x et PS=x*n
  65.          PS    = 0.D0
  66.          DO 210 ID = 1,IDIM
  67.             IDD1 = 3 + ID
  68.             XVALB(I,IND,IDD1) = XPTB(NPOI,1,ID)
  69.             PS = PS + XPTB(NPOI,1,ID) * XPALB(I,ID1+ID)
  70.  210     CONTINUE
  71. c        x" = deplacement dans le plan du cercle - excentration :
  72. c        x" = x'-e  avec  x'= x-(x*n)*n   et  XDEP = sqrt(x"*x")
  73. c        rem bp : excentration hors plan fausse l'enfoncement, non ?
  74.          PSXPME = 0.D0
  75.          DO 212 ID = 1,IDIM
  76.             XPRIM = XPTB(NPOI,1,ID) - PS * XPALB(I,ID1+ID)
  77.             XPME  = XPRIM - XPALB(I,ID2+ID)
  78.             PSXPME = PSXPME + XPME**2
  79.  212     CONTINUE
  80.          XDEP = SQRT(PSXPME)
  81.  
  82. *  calcul de la force de choc
  83. c       -cas sans amortissement
  84.          IF (ITYP.EQ.21) THEN
  85.             CALL DYCHE2(XDEP,XRAID,XJEU, XFL,iannul)
  86.  
  87. c       -cas avec amortissement : il faut la vitesse normale au choc
  88.          ELSE
  89.             XAMO   = XPALB(I,3)
  90. cbp,2020-09            XDEPM1 = 0.D0
  91.             XVIT = 0.D0
  92. *            IF (XDEP.GT.1D-20) THEN
  93.             IF (XDEP.GT.xpetit) THEN
  94.                PS2 = 0.D0
  95.                DO  213 ID=1,IDIM
  96.                   PS2 = PS2 + XPTB(NPOI,2,ID) * XPALB(I,ID1+ID)
  97.  213           CONTINUE
  98. cbp,2020-09 : intro de XVITN = vitesse hors plan du cercle (= selon n)
  99.                XVITN = SQRT(PS2)
  100.                DO 214 ID = 1,IDIM
  101. cbp,2020-09                  XPRIM = XPTB(NPOI,IND,ID) - PS * XPALB(I,ID1+ID)
  102. cbp,2020-09                  XPME = XPRIM - XPALB(I,ID2+ID)
  103. cbp,2020-09                  XNOR = XPME / XDEP
  104. cbp,2020-09                  XDEPM1 = XDEPM1 + (XPTB(NPOI,IND2,ID)
  105. cbp,2020-09     &             -PS2*XPALB(I,ID1+ID)-XPALB(I,ID2+ID))*XNOR
  106.                   XVIT = XVIT 
  107.      &            + (XPTB(NPOI,2,ID) - XVITN * XPALB(I,ID1+ID))**2
  108.  214           CONTINUE
  109.                XVIT = SQRT(XVIT)
  110.             ENDIF
  111. cbp,2020-09            XVIT = (XDEP - XDEPM1) / PDTS2
  112.             XVALB(I,IND,3) = XVIT
  113.             CALL DYCHA2(XDEP,XVIT,XRAID,XJEU,XAMO, XFL,iannul)
  114.          ENDIF
  115.  
  116. *  stockage
  117.          XVALB(I,IND,1) = XFL
  118. *         IF (XDEP.GE.XJEU .AND. XDEP.GT.1D-20) THEN
  119.          IF (XDEP.GE.XJEU.AND.XDEP.GT.xpetit) THEN
  120.             XPME = 0.D0
  121.             DO 216 ID = 1,IDIM
  122.                XPRIM = XPTB(NPOI,1,ID) - PS * XPALB(I,ID1+ID)
  123.                XPME = XPRIM - XPALB(I,ID2+ID)
  124.                XNOR = XPME / XDEP
  125.                FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFL * XNOR
  126.  216        CONTINUE
  127.          ENDIF
  128.  
  129.  
  130. *--------------------------------------------------------------------*
  131. * --- choc elementaire point_cercle_frottement
  132. *                      avec ou sans amortissement
  133. *                      interieur ou exterieur
  134. *--------------------------------------------------------------------*
  135. *
  136. cbp,2020       ELSE IF (ITYP.EQ.23  .OR. ITYP.EQ.24 
  137. cbp,2020      &    .or. ITYP.EQ.123 .OR. ITYP.EQ.124) THEN
  138.  
  139.       ELSE IF (ITYP.EQ.24 .OR. ITYP.EQ.124) THEN
  140.         write(*,*) 'BP: l amortissement ne change pas le numero de ITYP'
  141.         call erreur(5)
  142.         return
  143.  
  144.       ELSE IF (ITYP.EQ.23  .OR. ITYP.EQ.123 ) THEN
  145.  
  146.          NPOI   = IPLIB(I,1)
  147.          IGP    = IPALB(I,2)
  148.          IDIM   = IPALB(I,3)
  149.          if (ITYP.LT.100) then
  150.            INTER=1
  151.          else
  152.            INTER=0
  153.          endif
  154. cbp,2020          IF (ITYP.EQ.23 .or. ITYP.EQ.123) THEN
  155. cbp,2020             ID1 = 6
  156. cbp,2020          ELSE
  157. cbp,2020             ID1 = 7
  158. cbp,2020          ENDIF
  159.          ID1 = 10
  160.          ID2 = ID1 + IDIM
  161.          ID3 = ID1 + 2*IDIM
  162.          ID4 = ID1 + 3*IDIM
  163.          ID5 = ID1 + 4*IDIM
  164.          ID6 = ID1 + 5*IDIM
  165.          ID7 = ID1 + 6*IDIM
  166.          ID8 = ID1 + 7*IDIM
  167.          ID9 = ID1 + 8*IDIM
  168.  
  169. *        Vitesse d'entrainement tangentielle : Ve=Omega*R
  170.          Ve = XPALB(I,10)
  171. *  Attention : on considere un point (et pas un cercle) => on a : R=jeu!
  172. *              => comme R=jeu, les phases d'adherence sont celles d'un 
  173. *                 contact point_cercle et pas cercle_cercle
  174.          Omega = Ve / XPALB(I,2)
  175.  
  176. c        write(*,*) 'devfb3: ITYP,Ve,Omega,IGP,IND=',ITYP,Ve,Omega,IGP,IND
  177. *  si glissement au pas precedent, reactualisation de la position
  178. *  origine d'adherence
  179.          IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  180.             DO 230 ID=1,IDIM
  181.                XPALB(I,ID7+ID) = XPTB(NPOI,1,ID)
  182.  230        CONTINUE
  183. c        write(*,*) 'devfb3: x_adhe0=',(XPALB(I,ID7+iou),iou=1,IDIM),ID7
  184. cbp,2020 si adherence et entrainement, le point d'origine avance !
  185. c        (rem : on utilise une grande rotation)
  186.          ELSEIF (IND.ne.3) THEN
  187. c        write(*,*) 'devfb3: x_adhe0=',(XPALB(I,ID7+iou),iou=1,IDIM),ID7
  188. c           on se place en 3D uniquement (IDIM=3) 
  189. c           temporairement (!),
  190. c           ID4: direction du point d'adherence d'origine dans le plan du cercle
  191. c           ID5: direction orthogonale dans le plan du cercle
  192.             PS=0.D0
  193.             DO ID=1,3
  194.               PS=PS+XPALB(I,ID7+ID)*XPALB(I,ID1+ID)
  195.             ENDDO
  196.             PS4=0.D0
  197.             DO ID=1,3
  198.               XADH4=XPALB(I,ID7+ID)-PS*XPALB(I,ID1+ID)
  199.               XPALB(I,ID4+ID)=XADH4
  200.               PS4=PS4+XADH4*XADH4
  201.             ENDDO
  202. *         -cas ou le point d'adherence n'a pas ete defini
  203.            if(PS4.le.xpetit) then
  204.             DO 231 ID=1,IDIM
  205.                XPALB(I,ID7+ID) = XPTB(NPOI,1,ID)
  206.  231        CONTINUE
  207. *         -cas "normal"
  208.            else
  209.             Radhe4=SQRT(PS4)
  210.             DO ID=1,3
  211.               XPALB(I,ID4+ID)=XPALB(I,ID4+ID)/Radhe4
  212.             ENDDO
  213.             XPALB(I,ID5+1)=XPALB(I,ID1+2)*XPALB(I,ID4+3)
  214.      &                    -XPALB(I,ID1+3)*XPALB(I,ID4+2)
  215.             XPALB(I,ID5+2)=XPALB(I,ID1+3)*XPALB(I,ID4+1)
  216.      &                    -XPALB(I,ID1+1)*XPALB(I,ID4+3)
  217.             XPALB(I,ID5+3)=XPALB(I,ID1+1)*XPALB(I,ID4+2)
  218.      &                    -XPALB(I,ID1+2)*XPALB(I,ID4+1)
  219. c           rotation du point d'origine d'adherence
  220.             RCOS4=Radhe4*COS(Omega*PDTS2)
  221.             RSIN4=Radhe4*SIN(Omega*PDTS2)
  222.             DO ID=1,3
  223.                XPALB(I,ID7+ID)=RCOS4*XPALB(I,ID4+ID)
  224.      &                        +RSIN4*XPALB(I,ID5+ID)
  225.             ENDDO
  226.           endif
  227.          ENDIF
  228.  
  229. *  calcul du deplacement normal au plan du Cercle
  230. c        PS = x(t_ind) * nCercle
  231. *  et des variations de deplacement suivantes
  232. c        XPALB(I,ID4+ID) = | dx = x(t_ind)-x(t_ind2) (avant bp,2020-09)
  233. c                          | v                       (apres bp,2020-09) 
  234. c        XPALB(I,ID5+ID) = dx_adhe = x(t_ind)-x(t_adherence)
  235.          PS = 0.D0
  236.          DO 232 ID = 1,IDIM
  237.             IDD2 = 3 + ID
  238.             XVALB(I,IND,IDD2) = XPTB(NPOI,1,ID)
  239. cbp,2020-09              XPALB(I,ID4+ID) = XPTB(NPOI,IND,ID) - XPTB(NPOI,IND2,ID)
  240.             XPALB(I,ID4+ID) = XPTB(NPOI,2,ID)
  241.             XPALB(I,ID5+ID) = XPTB(NPOI,1,ID) - XPALB(I,ID7+ID)
  242.             PS = PS + XPTB(NPOI,1,ID) * XPALB(I,ID1+ID)
  243.  232     CONTINUE
  244.  
  245. *  calcul de la normale au plan tangent du contact, 
  246. c        XPRIM = x(t_ind) - (x(t_ind)*nCercle)*nCercle 
  247. c              = deplacement dans le plan du cercle
  248. c        XPALB(I,ID3+ID) = OX dans la boucle 234 
  249. c                        = AX (deplacement) - AO(excentrement)
  250. c                        avec O centre du cercle et X position du point
  251. c                   puis = n = normale de contact apres boucle 236
  252. *  et de la valeur du deplacement suivant cette normale
  253. c        XDEP = |OX|
  254.          XPME = 0.D0
  255.          PSXPME = 0.D0
  256.          DO 234 ID = 1,IDIM
  257.             XPRIM = XPTB(NPOI,1,ID) - PS * XPALB(I,ID1+ID)
  258.             XPME = XPRIM - XPALB(I,ID2+ID)
  259.             XPALB(I,ID3+ID) = XPME
  260.             PSXPME = PSXPME + XPME * XPME
  261.  234     CONTINUE
  262.          XDEP = SQRT(PSXPME)
  263. *         IF(XDEP.GT.1.D-20) THEN
  264.          IF(XDEP.GT.xpetit) THEN
  265.             DO 236 ID = 1,IDIM
  266.                XPALB(I,ID3+ID) = XPALB(I,ID3+ID) / XDEP
  267.  236        CONTINUE
  268.          ENDIF
  269.  
  270. *  calcul de la vitesse tangentielle par derivee a gauche
  271. c        XPALB(I,ID4+ID) =| (dx - (dx*n)*n)/dt (avant bp,2020-09)
  272. *                         |  v  -  (v*n)*n     (apres bp,2020-09)
  273. *  et la variation de deplacement par rapport a la position
  274. c  origine d'adherence sur le plan tangent
  275. c        XPALB(I,ID5+ID) = dx_adhe -  (dx_adhe*n)*n
  276.          XVITN= 0.D0
  277.          PSN0 = 0.D0
  278.          DO 238 ID = 1,IDIM
  279.             XVITN= XVITN+ XPALB(I,ID4+ID) * XPALB(I,ID3+ID)
  280.             PSN0 = PSN0 + XPALB(I,ID5+ID) * XPALB(I,ID3+ID)
  281.  238     CONTINUE
  282.          DO 240 ID = 1,IDIM
  283.             XPALB(I,ID4+ID) = XPALB(I,ID4+ID) - XVITN* XPALB(I,ID3+ID)
  284.             XPALB(I,ID5+ID) = XPALB(I,ID5+ID) - PSN0 * XPALB(I,ID3+ID)
  285.  240     CONTINUE
  286. c        vitesse normale et tangentielle
  287. cbp,2020-09            XVITN = XVITN / PDTS2
  288. cbp,2020-09            DO ID = 1,IDIM
  289. cbp,2020-09               XPALB(I,ID4+ID) = XPALB(I,ID4+ID) / PDTS2
  290. cbp,2020-09            ENDDO
  291. cbp,2020 calcul de la vitesse d'entrainement = Ve*(nCercle pvec n)
  292. c        temporairement, ID6: direction de la vitesse d'entrainement 
  293.             XPALB(I,ID6+1)=XPALB(I,ID1+2)*XPALB(I,ID3+3)
  294.      &                    -XPALB(I,ID1+3)*XPALB(I,ID3+2)
  295.             XPALB(I,ID6+2)=XPALB(I,ID1+3)*XPALB(I,ID3+1)
  296.      &                    -XPALB(I,ID1+1)*XPALB(I,ID3+3)
  297.             XPALB(I,ID6+3)=XPALB(I,ID1+1)*XPALB(I,ID3+2)
  298.      &                    -XPALB(I,ID1+2)*XPALB(I,ID3+1)
  299. cbp,2020 vitesse relative = absolue - entrainement  
  300.          DO ID=1,IDIM
  301.             XPALB(I,ID4+ID)=XPALB(I,ID4+ID) - Ve*XPALB(I,ID6+ID)
  302.          ENDDO
  303. c          write(*,*) 'devfb3: NORMALE=',(XPALB(I,ID1+iou),iou=1,IDIM),ID1
  304. c          write(*,*) 'devfb3: EXCENTR=',(XPALB(I,ID2+iou),iou=1,IDIM),ID2
  305. c          write(*,*) 'devfb3: n_choc=',(XPALB(I,ID3+iou),iou=1,IDIM),ID3
  306. c          write(*,*) 'devfb3: Vt_rel=',(XPALB(I,ID4+iou),iou=1,IDIM),ID4
  307. c          write(*,*) 'devfb3: dx_adh=',(XPALB(I,ID5+iou),iou=1,IDIM),ID5
  308. cbp,2020 ID6 reprendra son role de Ft dans DYCHA31
  309.  
  310. *  calcul de la force de choc
  311. cbp,2020c       -cas sans amortissement
  312. cbp,2020         IF (ITYP.EQ.23 .or. ITYP.EQ.123) THEN
  313. cbp,2020            CALL DYCHE3(XDEP,IDIM,IGP,XPALB,NLIAB,I,INTER,
  314. cbp,2020     &                  XFN,XFT,XPUS,iannul)
  315.  
  316. c       -cas avec amortissement : il faut la vitesse normale
  317. cbp,2020         ELSE
  318. cbp,2020-09        XVITN = PSN / PDTS2
  319.             XVALB(I,IND,3) = XVITN
  320. cbp,2020            CALL DYCHA3(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER
  321.             CALL DYCHA31(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER
  322.      &                  ,XFN,XFT,XPUS,iannul)
  323. cbp,2020         ENDIF
  324.          XVALB(I,IND,1) = XFN
  325.          XVALB(I,IND,10) = ABS(XFT)
  326.          XVALB(I,IND,12) = XPUS
  327.          IPALB(I,2) = IGP
  328. c          write(*,*) 'devfb3: Xn,Vn,Fn=',XDEP,XVITN,XFN
  329.  
  330. *  stockage
  331.          IF (IGP .EQ. 1)  THEN
  332.             PS = 0.D0
  333.               DO 20 ID = 1,IDIM
  334.                 PS = PS + (XPALB(I,ID4+ID)*XPALB(I,ID4+ID))
  335.  20           CONTINUE
  336.             XVITT = SQRT(PS)
  337.          ELSE
  338.             XVITT = 0.D0
  339.          ENDIF
  340.          XVALB(I,IND,11) = XVITT
  341. *  si glissement, memorisation de la vitesse tangentielle et de la force
  342. *  tangentielle
  343.          IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  344.             DO 242 ID = 1,IDIM
  345.                XPALB(I,ID8+ID) = XPALB(I,ID4+ID)
  346.                XPALB(I,ID9+ID) = XPALB(I,ID6+ID)
  347.  242        CONTINUE
  348.          ENDIF
  349.          DO 244 ID = 1,IDIM
  350.             FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFN*XPALB(I,ID3+ID)
  351.      &                                      + XPALB(I,ID6+ID)
  352.  244     CONTINUE
  353.  
  354.  
  355. *--------------------------------------------------------------------*
  356. * --- choc elementaire point_cercle_mobile
  357. *                      avec ou sans amortissement
  358. *--------------------------------------------------------------------*
  359. * on neglige la rotation (torsion) du cercle
  360.  
  361. cbp      ELSE IF (ITYP.EQ.33 .OR. ITYP.EQ.34) THEN
  362.       ELSE IF (ITYP.EQ.33  .OR. ITYP.EQ.34 
  363.      &    .or. ITYP.EQ.133 .OR. ITYP.EQ.134) THEN
  364.          NPOA   = IPLIB(I,1)
  365.          NPOB   = IPLIB(I,2)
  366.          IGP    = IPALB(I,2)
  367.          IDIM   = IPALB(I,3)
  368. cbp         INTER  = IPALB(I,4)
  369.          if (ITYP.LT.100) then
  370.            INTER=1
  371.          else
  372.            INTER=0
  373.          endif
  374.          IF (ITYP.EQ.33 .or. ITYP.EQ.133) THEN
  375.             ID1 = 6
  376.          ELSE
  377.             ID1 = 7
  378.          ENDIF
  379.          ID2 = ID1 + IDIM
  380.          ID3 = ID1 + 2*IDIM
  381.          ID4 = ID1 + 3*IDIM
  382.          ID5 = ID1 + 4*IDIM
  383.          ID6 = ID1 + 5*IDIM
  384.          ID7 = ID1 + 6*IDIM
  385.          ID8 = ID1 + 7*IDIM
  386.          ID9 = ID1 + 8*IDIM
  387. *  si pas d'adherence (c.a.d. glissement ou pas de contact) au pas precedent,
  388. *  reactualisation de la position ecart origine d'adherence
  389.          IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  390.             DO 330 ID=1,IDIM
  391.                XPALB(I,ID7+ID) = (XPTB(NPOA,1,ID) +
  392.      &                           FEXPSM(NPOA,NPA,IND1,ID) )
  393.      &                       - ( XPTB(NPOB,1,ID) +
  394.      &                           FEXPSM(NPOB,NPA,IND1,ID))
  395.  330        CONTINUE
  396.          ENDIF
  397. *  calcul du deplacement courant relatif (PS) sur la normale au cercle
  398.          PS = 0.D0
  399.          DO 332 ID = 1,IDIM
  400.             IDD2 = 3 + ID
  401. c           x_n = xA_n - xB_n
  402.             XDE2 = XPTB(NPOA,1,ID) - XPTB(NPOB,1,ID)
  403.             XDE2 = XDE2 + FEXPSM(NPOa,NPA,IND1,ID)
  404.      &                  - FEXPSM(NPOb,NPA,IND1,ID)
  405.             XVALB(I,IND,IDD2) = XDE2
  406. cbp,2020-09c           x_n-1/2 = xA_n-1/2 - xB_n-1/2
  407. cbp,2020-09            XDM2 = XPTB(NPOA,IND2,ID) - xptb (npob,ind2,id)
  408. cbp,2020-09            XDM2 = XDM2 + FEXPSM(NPOA,NPAM1,INDM1,ID)
  409. cbp,2020-09     &                  - FEXPSM(NPOb,NPAM1,INDM1,ID)
  410. cbp,2020-09c           dx = x_n - x_n-1/2
  411. cbp,2020-09            XPALB(I,ID4+ID) = XDE2 - XDM2
  412. c           v_n = vA_n - vB_n
  413.             XPALB(I,ID4+ID) = XPTB(NPOA,2,ID) - XPTB(NPOB,2,ID)
  414. c           dx_adhe = x_n - x_adherence0
  415.             XPALB(I,ID5+ID) = XDE2 - XPALB(I,ID7+ID)
  416. c           PS = x_n * nCercle 
  417.             PS = PS + XDE2 * XPALB(I,ID1+ID)
  418.  332     CONTINUE
  419.  
  420. *  calcul de la valeur du deplacement relatif suivant la normale
  421. *  au plan tangent du contact (on tient compte de l'excentrement)
  422.          XPME = 0.D0
  423.          PSXPME = 0.D0
  424.          DO 334 ID = 1,IDIM
  425. c           x_n = xA_n - xB_n
  426.             XDE2 = XPTB(NPOA,1,ID) - XPTB(NPOB,1,ID)
  427.             XDE2 = XDE2 + FEXPSM(NPOa,NPA,IND1,ID)
  428.      &                  - FEXPSM(NPOb,NPA,IND1,ID)
  429. c           x' = x_n - (x_n * nCercle) * nCercle
  430.             XPRIM = XDE2 - PS * XPALB(I,ID1+ID)
  431. c           x" = x'-e  et  XDEP = sqrt(x"*x")
  432. c           rem bp : excentration hors plan fausse l'enfoncement, non ?
  433.             XPME = XPRIM - (XPALB(I,ID2+ID) -
  434.      &                     XPALB(I,ID2+ID)*XPALB(I,ID1+ID))
  435.             XPALB(I,ID3+ID) = XPME
  436.             PSXPME = PSXPME + XPME**2
  437.  334     CONTINUE
  438.          XDEP = SQRT(PSXPME)
  439.  
  440. * normale au plan tangent = direction normale du choc = n
  441. *         IF(XDEP.GT.1.D-20) THEN
  442.          IF(XDEP.GT.xpetit) THEN
  443.             DO 336 ID = 1,IDIM
  444.                XPALB(I,ID3+ID) = XPALB(I,ID3+ID) / XDEP
  445.  336        CONTINUE
  446.          ENDIF
  447. *  calcul de la vitesse relative tangentielle : vt = v - (v*n)*n
  448. *  et de l'allongement d'adherence
  449.          XVITN = 0.D0
  450.          PSN0  = 0.D0
  451.          DO 338 ID = 1,IDIM
  452.             XVITN= XVITN+ XPALB(I,ID4+ID) * XPALB(I,ID3+ID)
  453.             PSN0 = PSN0 + XPALB(I,ID5+ID) * XPALB(I,ID3+ID)
  454.  338     CONTINUE
  455.          DO 340 ID = 1,IDIM
  456.             XPALB(I,ID4+ID) = (XPALB(I,ID4+ID) - XVITN* XPALB(I,ID3+ID))
  457. cbp,2020-09     &                         /PDTS2
  458.             XPALB(I,ID5+ID) =  XPALB(I,ID5+ID) - PSN0 * XPALB(I,ID3+ID)
  459.  340     CONTINUE
  460. *  calcul de la force de choc
  461.          IF (ITYP.EQ.33 .or. ITYP.EQ.133) THEN
  462.             CALL DYCHE3(XDEP,IDIM,IGP,XPALB,NLIAB,I,INTER
  463.      &                 ,XFN,XFT,XPUS,iannul)
  464.          ELSE
  465. cbp,2020-09            XVITN = PSN / PDTS2
  466.             XVALB(I,IND,3) = XVITN
  467.             CALL DYCHA3(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER
  468.      &                 ,XFN,XFT,XPUS,iannul)
  469.          ENDIF
  470.          XVALB(I,IND,1) = XFN
  471.          XVALB(I,IND,10) = ABS(XFT)
  472.          XVALB(I,IND,12) = XPUS
  473.          IPALB(I,2) = IGP
  474.  
  475.          IF (IGP .EQ. 1)  THEN
  476.             PS = 0.D0
  477.               DO 30 ID = 1,IDIM
  478.                 PS = PS + (XPALB(I,ID4+ID)*XPALB(I,ID4+ID))
  479.  30           CONTINUE
  480.             XVITT = SQRT(PS)
  481.          ELSE
  482.             XVITT = 0.D0
  483.          ENDIF
  484.          XVALB(I,IND,11) = XVITT
  485. *  si pas d'adherence (c.a.d. glissement ou pas de contact)
  486. *  memorisation de la vitesse tangentielle et de la force tangentielle
  487.          IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  488.             DO 342 ID = 1,IDIM
  489.                XPALB(I,ID8+ID) = XPALB(I,ID4+ID)
  490.                XPALB(I,ID9+ID) = XPALB(I,ID6+ID)
  491.  342        CONTINUE
  492.          ENDIF
  493. *        fA = fA + |fn|*n + ft
  494. *        fB = fB - |fn|*n + ft
  495.          DO 344 ID = 1,IDIM
  496.             FTOTB(NPOA,ID) = FTOTB(NPOA,ID) + XFN*XPALB(I,ID3+ID)
  497.      &                                      + XPALB(I,ID6+ID)
  498.             FTOTB(NPOB,ID) = FTOTB(NPOB,ID) - XFN*XPALB(I,ID3+ID)
  499.      &                                      - XPALB(I,ID6+ID)
  500.  344     CONTINUE
  501.  
  502.  
  503. *--------------------------------------------------------------------*
  504. * --- choc elementaire cercle_cercle_frottement
  505. *                      avec ou sans amortissement
  506. *--------------------------------------------------------------------*
  507. *
  508. cbp      ELSE IF (ITYP.EQ.25 .OR. ITYP.EQ.26) THEN
  509.       ELSE IF (ITYP.EQ.25  .OR. ITYP.EQ.26 
  510.      &    .or. ITYP.EQ.125 .OR. ITYP.EQ.126) THEN
  511.          NPOI   = IPLIB(I,1)
  512.          IGP    = IPALB(I,2)
  513.          IDIM   = IPALB(I,3)
  514. cbp         INTER  = IPALB(I,4)
  515.          if (ITYP.LT.100) then
  516.            INTER=1
  517.          else
  518.            INTER=0
  519.          endif
  520.          IF (ITYP.EQ.25 .or. ITYP.EQ.125) THEN
  521.             ID1  = 6
  522.          ELSE
  523.             ID1  = 7
  524.          ENDIF
  525.          ID2  = ID1 + IDIM
  526.          ID3  = ID1 + 2*IDIM
  527.          ID4  = ID1 + 3*IDIM
  528.          ID5  = ID1 + 4*IDIM
  529.          ID6  = ID1 + 5*IDIM
  530.          ID7  = ID1 + 6*IDIM
  531.          ID8  = ID1 + 7*IDIM
  532.          ID9  = ID1 + 8*IDIM
  533.          ID10 = ID1 + 9*IDIM
  534. *        XRAYT : rayon du tube interne (rTube)= RAYON_SUPPORT
  535. *        XREXT : rayon du tube externe (rExt) = RAYON_BUTEE
  536.          XRAYT = XPALB(I,ID10+1)
  537.          XREXT = XPALB(I,2)
  538. *  calcul du deplacement du point fibre neutre dans le plan du cercle
  539. *  recuperation de la normale de choc au pas precedent
  540.          DO 249 ID = 1,IDIM
  541.             XN2(ID)= XPALB(I,ID3+ID)
  542.  249     CONTINUE
  543. *  calcul de la normale de choc au pas courant : XPALB(I,ID3+ID) = n
  544. *        PSXPN = x * nCercle
  545.          PSXPN = 0.D0
  546.          DO 250 ID = 1,IDIM
  547.             PSXPN = PSXPN + ( XPTB(NPOI,1,ID) * XPALB(I,ID1+ID) )
  548.  250     CONTINUE
  549.          PSXPME = 0.D0
  550.          DO 254 ID = 1,IDIM
  551. *           x" = x - (x*nCercle)*nCercle - exce
  552.             XXPME = ( XPTB(NPOI,1,ID) - ( PSXPN * XPALB(I,ID1+ID) ) )
  553.      &                                  - XPALB(I,ID2+ID)
  554.             XPALB(I,ID3+ID) = XXPME
  555.             PSXPME = PSXPME + XXPME**2
  556.  254     CONTINUE
  557. *        deplacement radial : PSXPME = |x"| = sqrt(x"*x")
  558. *        normale de choc    : XPALB(I,ID3+ID) = n = x"/sqrt(x"*x")
  559.          PSXPME = SQRT(PSXPME)
  560. *         IF (PSXPME.GT.1D-20) THEN
  561.          IF (PSXPME.GT.xpetit) THEN
  562.             DO 256 ID = 1,IDIM
  563.                XPALB(I,ID3+ID) = XPALB(I,ID3+ID) / PSXPME
  564.  256        CONTINUE
  565.          ENDIF
  566. *        deplacement du point geometrique de contact suivant la normale n 
  567. *        |xP| = |x"| + R
  568.          XDEP =  PSXPME + XRAYT
  569. *        deplacement du point geometrique de contact 
  570. *        xP = x + R*n
  571.          XPTP2(1) =  XPTB(NPOI,1,1) + (XPALB(I,ID3+1)*XRAYT)
  572.          XPTP2(2) =  XPTB(NPOI,1,2) + (XPALB(I,ID3+2)*XRAYT)
  573.          XPTP2(3) =  XPTB(NPOI,1,3) + (XPALB(I,ID3+3)*XRAYT)
  574. *  calcul du deplacement du point de contact au pas precedent
  575. cbp,2020-09         XPTPM2(1) =  XPTB(NPOI,IND2,1) + (XN2(1)*XRAYT)
  576. cbp,2020-09         XPTPM2(2) =  XPTB(NPOI,IND2,2) + (XN2(2)*XRAYT)
  577. cbp,2020-09         XPTPM2(3) =  XPTB(NPOI,IND2,3) + (XN2(3)*XRAYT)
  578. cbp,2020-09 : pour l'instant, on fait le calcul avec formule ci-dessous : 
  579. c a modifier + tard si schema integration change la relation deplacement<->vitesse
  580.          XPTPM2(1) = XPTB(NPOI,1,1)-PDTS2*XPTB(NPOI,2,1) 
  581.      &             + (XN2(1)*XRAYT)
  582.          XPTPM2(2) = XPTB(NPOI,1,1)-PDTS2*XPTB(NPOI,2,2) 
  583.      &             + (XN2(2)*XRAYT)
  584.          XPTPM2(3) = XPTB(NPOI,1,1)-PDTS2*XPTB(NPOI,2,3) 
  585.      &             + (XN2(3)*XRAYT)
  586. *  Vitesse *(-1) du point de contact appartenant a la structure mobile
  587. *  due a la rotation absolue
  588. cbp,2020-09               XVPC0(1)        =  (1.D0/ PDTS2) *
  589. cbp,2020-09     &             (  ( XPTB(NPOI,IND2,5) * XPALB(I,ID3+3) * XRAYT ) -
  590. cbp,2020-09     &                ( XPTB(NPOI,IND2,6) * XPALB(I,ID3+2) * XRAYT ) -
  591. cbp,2020-09     &                ( XPTB(NPOI,IND ,5) * XPALB(I,ID3+3) * XRAYT ) +
  592. cbp,2020-09     &                ( XPTB(NPOI,IND ,6) * XPALB(I,ID3+2) * XRAYT ) )
  593. cbp,2020-09               XVPC0(2)        =  (1.D0/ PDTS2) *
  594. cbp,2020-09     &             (  ( XPTB(NPOI,IND2,6) * XPALB(I,ID3+1) * XRAYT ) -
  595. cbp,2020-09     &                ( XPTB(NPOI,IND2,4) * XPALB(I,ID3+3) * XRAYT ) -
  596. cbp,2020-09     &                ( XPTB(NPOI,IND ,6) * XPALB(I,ID3+1) * XRAYT ) +
  597. cbp,2020-09     &                ( XPTB(NPOI,IND ,4) * XPALB(I,ID3+3) * XRAYT ) )
  598. cbp,2020-09               XVPC0(3)        =  (1.D0/ PDTS2) *
  599. cbp,2020-09     &             (  ( XPTB(NPOI,IND2,4) * XPALB(I,ID3+2) * XRAYT ) -
  600. cbp,2020-09     &                ( XPTB(NPOI,IND2,5) * XPALB(I,ID3+1) * XRAYT ) -
  601. cbp,2020-09     &                ( XPTB(NPOI,IND ,4) * XPALB(I,ID3+2) * XRAYT ) +
  602. cbp,2020-09     &                ( XPTB(NPOI,IND ,5) * XPALB(I,ID3+1) * XRAYT ) )
  603. *         -Vrota  = -w  x (R*n)      (x designe le produit vectoriel ∧)
  604.           XVPC0(1)= - (XPTB(NPOI,1 ,5) * XRAYT*XPALB(I,ID3+3))
  605.      &              + (XPTB(NPOI,1 ,6) * XRAYT*XPALB(I,ID3+2))
  606.           XVPC0(2)= - (XPTB(NPOI,1 ,6) * XRAYT*XPALB(I,ID3+1))
  607.      &              + (XPTB(NPOI,1 ,4) * XRAYT*XPALB(I,ID3+3))
  608.           XVPC0(3)= - (XPTB(NPOI,1 ,4) * XRAYT*XPALB(I,ID3+2))
  609.      &              + (XPTB(NPOI,1 ,5) * XRAYT*XPALB(I,ID3+1))
  610. *        si glissement au pas precedent, reactualisation de la position
  611. *        origine d'adherence a l'aide du point de contact
  612.          IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  613.             DO 257 ID=1,IDIM
  614.                XPALB(I,ID7+ID) = XPTP2(ID)
  615.  257        CONTINUE
  616. *        si adherence au pas precedent, on met a jour x(t_adhe0)
  617.          ELSE
  618. c           composante hors plan (selon nCercle) de Vrota
  619.             PXVPC0 = 0.D0
  620.             DO 258 ID = 1,IDIM
  621.                 PXVPC0 = PXVPC0 + XVPC0(ID)*XPALB(I,ID1+ID)
  622.  258        CONTINUE
  623. c           XVPCT : Vrota' = Vrota - (Vrota*nCercle)*nCercle
  624.             DO 259 ID = 1,IDIM
  625.                 XVPCT(ID) = XVPC0(ID) - PXVPC0 * XPALB(I,ID1+ID)
  626.  259        CONTINUE
  627. cbp,2020-09 : formule ci-dessous non comprise...???
  628.             XPALB(I,ID7+1) = XPALB(I,ID7+1)  +
  629.      &                  (XVPCT(1)*PDTS2)*(XREXT/(XREXT-XRAYT))
  630.             XPALB(I,ID7+2) = XPALB(I,ID7+2)  +
  631.      &                  (XVPCT(2)*PDTS2)*(XREXT/(XREXT-XRAYT))
  632.             XPALB(I,ID7+3) = XPALB(I,ID7+3)  +
  633.      &                  (XVPCT(3)*PDTS2)*(XREXT/(XREXT-XRAYT))
  634.          ENDIF
  635. *  calcul du deplacement sur la normale au plan de section droite
  636. *  et de l'ecart a la position origine adherencee
  637.          DO 260 ID = 1,IDIM
  638.             IDD1 =  3 + ID
  639.             IDD2 =  6 + ID
  640.             IDD3 = 15 + ID
  641. c           x : ddl de TRANSLATION au point A
  642.             XVALB(I,IND,IDD1) = XPTB(NPOI,1,ID)
  643. c           vP : vitesse du point P de contact (a modifier + tard)
  644.             XVALB(I,IND,IDD2) = (XPTP2(ID) - XPTPM2(ID) ) / PDTS2
  645. c           theta : ddl de ROTATION au point A
  646.             XVALB(I,IND,IDD3) = XPTB(NPOI,1,ID+3)
  647. c           mvt du point geometrique de contact - mvt du point de contact attache au tube
  648. c           in fine, on a vt = dxP/dt + Vrota
  649. cbp,2020-09 XPALB(I,ID4+ID) = XPTP2(ID) - XPTPM2(ID) - (XVPC0(ID)*PDTS2)
  650.             XPALB(I,ID4+ID) = (XPTP2(ID)-XPTPM2(ID))/PDTS2 - XVPC0(ID)
  651. c           allongement du ressort d'adherence (pas compris...??? 
  652. c           notamment, pourquoi different du glissement?)
  653.             XPALB(I,ID5+ID) = XPTP2(ID) - XPALB(I,ID7+ID)
  654.  260     CONTINUE
  655. *  calcul de la vitesse tangentielle 
  656. *  et de l'ecart a la position origine adherencee dans le plan
  657.          XVITN= 0.D0
  658.          PSN0 = 0.D0
  659.          DO 262 ID = 1,IDIM
  660.             XVITN= XVITN+ XPALB(I,ID4+ID) * XPALB(I,ID3+ID)
  661.             PSN0 = PSN0 + XPALB(I,ID5+ID) * XPALB(I,ID3+ID)
  662.  262     CONTINUE
  663.          DO 264 ID = 1,IDIM
  664.             XPALB(I,ID4+ID) = XPALB(I,ID4+ID) - XVITN * XPALB(I,ID3+ID)
  665. cbp,2020-09     &                        (...)/PDTS2
  666.             XPALB(I,ID5+ID) = XPALB(I,ID5+ID) - PSN0  * XPALB(I,ID3+ID)
  667.  264     CONTINUE
  668. *  calcul de la force de choc
  669.          IF (ITYP.EQ.25 .or. ITYP.EQ.125) THEN
  670.             CALL DYCHE3(XDEP,IDIM,IGP,XPALB,NLIAB,I,INTER
  671.      &                  ,XFN,XFT,XPUS,iannul)
  672.          ELSE
  673.             XVITN = PSN / PDTS2
  674.             XVALB(I,IND,3) = XVITN
  675.             CALL DYCHA3(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER
  676.      &                  ,XFN,XFT,XPUS,iannul)
  677.          ENDIF
  678.          XVALB(I,IND,1) = XFN
  679.          XVALB(I,IND,10) = ABS(XFT)
  680.          XVALB(I,IND,12) = XPUS
  681.          IPALB(I,2) = IGP
  682. *  si glissement, memorisation de la vitesse tangentielle et de la force
  683. *  tangentielle
  684.          IF (IGP.EQ.1) THEN
  685.             DO 266 ID = 1,IDIM
  686.                XPALB(I,ID8+ID) = XPALB(I,ID4+ID)
  687.                XPALB(I,ID9+ID) = XPALB(I,ID6+ID)
  688.  266           CONTINUE
  689.          ELSE
  690.             DO 267 ID = 1,IDIM
  691.                XPALB(I,ID9+ID) = 0.D0
  692.  267        CONTINUE
  693.          ENDIF
  694. *        Force : f = fn*n + ft
  695.          DO 268 ID = 1,IDIM
  696.             XFOR = ( XFN * XPALB(I,ID3+ID) ) + XPALB(I,ID6+ID)
  697.             FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFOR
  698.             XFNT(ID) = XPALB (I ,ID6+ID)
  699.  268     CONTINUE
  700. *        Moment : m = (R*n) x f
  701.          XAPP1 = XRAYT * XPALB(I,ID3+1)
  702.          XAPP2 = XRAYT * XPALB(I,ID3+2)
  703.          XAPP3 = XRAYT * XPALB(I,ID3+3)
  704.          XAPFP1 = ( XAPP2 * XFNT(3) ) - ( XAPP3 * XFNT(2) )
  705.          XAPFP2 = ( XAPP3 * XFNT(1) ) - ( XAPP1 * XFNT(3) )
  706.          XAPFP3 = ( XAPP1 * XFNT(2) ) - ( XAPP2 * XFNT(1) )
  707.          XVALB(I,IND,13) = XAPFP1
  708.          XVALB(I,IND,14) = XAPFP2
  709.          XVALB(I,IND,15) = XAPFP3
  710.          FTOTB(NPOI,4) = FTOTB(NPOI,4) + XAPFP1
  711.          FTOTB(NPOI,5) = FTOTB(NPOI,5) + XAPFP2
  712.          FTOTB(NPOI,6) = FTOTB(NPOI,6) + XAPFP3
  713.  
  714.  
  715. *--------------------------------------------------------------------*
  716. * --- choc ...........
  717. *--------------------------------------------------------------------*
  718. *
  719. *     else if (ityp.eq.  ) then
  720. *        .......
  721. *        .......
  722. *
  723.       ENDIF
  724. *
  725.       END
  726.  
  727.  
  728.  
  729.  
  730.  
  731.  
  732.  
  733.  
  734.  
  735.  
  736.  
  737.  

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