C DEVFB3 SOURCE CHAT 05/01/12 22:45:25 5004 SUBROUTINE DEVFB3(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB, & NPLB,IND,IND1,INDM1,NPA,NPAM1,IND2,PDT,PDTS2, & FEXPSM,NPC1,I,iannul) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) *--------------------------------------------------------------------* * * * operateur dyne : algorithme de fu - de vogelaere * * ________________________________________________ * * * * calcul des forces de choc sur base b pour les liaisons * * de type point_cercle. * * * * param}tres: * * * * e ityp type de la liaison. * * es ftotb forces exterieures totalisees sur la base b. * * e xptb tableau des deplacements des points * * e ipalb renseigne sur la liaison. * * e iplib tableau contenant les numeros "dyne" de la liaison. * * e xpalb tableau contenant les parametres de la liaison. * * es xvalb tableau contenant les variables internes de liaisons. * * e nliab nombre de liaisons sur la base b. * * e nplb nombre total de points intervenant dans les liaisons. * * e ind indice du pas. * * e i num{ro de la liaison. * * * * * * auteur, date de creation: * * * * lionel vivan : le 22 septembre 1989 : creation * * bertrand beaufils : le 31 mai 1990 : ajout frottement sec* * I. Politopoulos : juin 1996 : ajout point_cercle_mobile * * * *--------------------------------------------------------------------* * -INC CCREEL INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*) REAL*8 XPALB(NLIAB,*),XPTB(NPLB,4,*),FTOTB(NPLB,*) REAL*8 XVALB(NLIAB,4,*),FEXPSM(NPLB,NPC1,2,*) REAL*8 XPTP2(3),XPTPM2(3),XFNT(3) REAL*8 XVPC0(3),XN2(3),XVPCT(3) * * --- choc elementaire point_cercle * avec ou sans amortissement * IF (ITYP.EQ.21 .OR. ITYP.EQ.22) THEN NPOI = IPLIB(I,1) IDIM = IPALB(I,3) XRAID = XPALB(I,1) XJEU = XPALB(I,2) IF (ITYP.EQ.21) THEN ID1 = 2 ELSE XAMO = XPALB(I,3) ID1 = 3 ENDIF ID2 = ID1 + IDIM PS = 0.D0 DO 210 ID = 1,IDIM IDD1 = 3 + ID XVALB(I,IND,IDD1) = XPTB(NPOI,IND,ID) PS = PS + XPTB(NPOI,IND,ID) * XPALB(I,ID1+ID) 210 CONTINUE * end do XPME = 0.D0 PSXPME = 0.D0 DO 212 ID = 1,IDIM XPRIM = XPTB(NPOI,IND,ID) - PS * XPALB(I,ID1+ID) XPME = XPRIM - XPALB(I,ID2+ID) PSXPME = PSXPME + XPME * XPME 212 CONTINUE * end do XDEP = SQRT(PSXPME) IF (ITYP.EQ.21) THEN CALL DYCHE2(XDEP,XRAID,XJEU, XFL,iannul) ELSE XDEPM1 = 0.D0 * IF (XDEP.GT.1D-20) THEN IF (XDEP.GT.xpetit) THEN PS2 = 0.D0 * correction gayffier 4/04/96 xdepm1 mal calculé DO 213 ID=1,IDIM PS2 = PS2 + XPTB(NPOI,IND2,ID) * XPALB(I,ID1+ID) 213 CONTINUE DO 214 ID = 1,IDIM XPRIM = XPTB(NPOI,IND,ID) - PS * XPALB(I,ID1+ID) XPME = XPRIM - XPALB(I,ID2+ID) XNOR = XPME / XDEP XDEPM1 = XDEPM1 + (XPTB(NPOI,IND2,ID) & -PS2*XPALB(I,ID1+ID)-XPALB(I,ID2+ID))*XNOR 214 CONTINUE * end do ENDIF XVIT = (XDEP - XDEPM1) / PDTS2 XVALB(I,IND,3) = XVIT CALL DYCHA2(XDEP,XVIT,XRAID,XJEU,XAMO, XFL,iannul) ENDIF XVALB(I,IND,1) = XFL * IF (XDEP.GE.XJEU .AND. XDEP.GT.1D-20) THEN IF (XDEP.GE.XJEU.AND.XDEP.GT.xpetit) THEN XPME = 0.D0 DO 216 ID = 1,IDIM XPRIM = XPTB(NPOI,IND,ID) - PS * XPALB(I,ID1+ID) XPME = XPRIM - XPALB(I,ID2+ID) XNOR = XPME / XDEP FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFL * XNOR 216 CONTINUE * end do ENDIF * * --- choc elementaire point_cercle_frottement * avec ou sans amortissement * ELSE IF (ITYP.EQ.23 .OR. ITYP.EQ.24) THEN NPOI = IPLIB(I,1) IGP = IPALB(I,2) IDIM = IPALB(I,3) INTER = IPALB(I,4) IF (ITYP.EQ.23) THEN ID1 = 6 ELSE ID1 = 7 ENDIF ID2 = ID1 + IDIM ID3 = ID1 + 2*IDIM ID4 = ID1 + 3*IDIM ID5 = ID1 + 4*IDIM ID6 = ID1 + 5*IDIM ID7 = ID1 + 6*IDIM ID8 = ID1 + 7*IDIM ID9 = ID1 + 8*IDIM * si glissement au pas precedent, reactualisation de la position * origine d'adherence IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 230 ID=1,IDIM XPALB(I,ID7+ID) = XPTB(NPOI,IND,ID) 230 CONTINUE * end do ENDIF * calcul du deplacement sur la normale au plan de section droite PS = 0.D0 DO 232 ID = 1,IDIM IDD2 = 3 + ID XVALB(I,IND,IDD2) = XPTB(NPOI,IND,ID) XPALB(I,ID4+ID) = XPTB(NPOI,IND,ID) - XPTB(NPOI,IND2,ID) XPALB(I,ID5+ID) = XPTB(NPOI,IND,ID) - XPALB(I,ID7+ID) PS = PS + XPTB(NPOI,IND,ID) * XPALB(I,ID1+ID) 232 CONTINUE * end do * calcul de la normale au plan tangent du contact, et de la valeur du * deplacement suivant cette normale XPME = 0.D0 PSXPME = 0.D0 DO 234 ID = 1,IDIM XPRIM = XPTB(NPOI,IND,ID) - PS * XPALB(I,ID1+ID) XPME = XPRIM - XPALB(I,ID2+ID) XPALB(I,ID3+ID) = XPME PSXPME = PSXPME + XPME * XPME 234 CONTINUE * end do XDEP = SQRT(PSXPME) * IF(XDEP.GT.1.D-20) THEN IF(XDEP.GT.xpetit) THEN DO 236 ID = 1,IDIM XPALB(I,ID3+ID) = XPALB(I,ID3+ID) / XDEP 236 CONTINUE * end do ENDIF * calcul de la vitesse tangentielle par derivee a gauche PSN = 0.D0 PSN0 = 0.D0 DO 238 ID = 1,IDIM PSN = PSN + XPALB(I,ID4+ID) * XPALB(I,ID3+ID) PSN0 = PSN0 + XPALB(I,ID5+ID) * XPALB(I,ID3+ID) 238 CONTINUE * end do DO 240 ID = 1,IDIM XPALB(I,ID4+ID) = (XPALB(I,ID4+ID) - & PSN * XPALB(I,ID3+ID))/PDTS2 XPALB(I,ID5+ID) = XPALB(I,ID5+ID) - & PSN0 * XPALB(I,ID3+ID) 240 CONTINUE * end do * calcul de la force de choc IF (ITYP.EQ.23) THEN CALL DYCHE3(XDEP,IDIM,IGP,XPALB,NLIAB,I,INTER, & XFN,XFT,XPUS,iannul) ELSE XVITN = PSN / PDTS2 XVALB(I,IND,3) = XVITN CALL DYCHA3(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER & ,XFN,XFT,XPUS,iannul) ENDIF XVALB(I,IND,1) = XFN XVALB(I,IND,10) = ABS(XFT) XVALB(I,IND,12) = XPUS IPALB(I,2) = IGP IF (IGP .EQ. 1) THEN PS = 0.D0 DO 20 ID = 1,IDIM PS = PS + (XPALB(I,ID4+ID)*XPALB(I,ID4+ID)) 20 CONTINUE XVITT = SQRT(PS) ELSE XVITT = 0.D0 ENDIF XVALB(I,IND,11) = XVITT * si glissement, memorisation de la vitesse tangentielle et de la force * tangentielle IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 242 ID = 1,IDIM XPALB(I,ID8+ID) = XPALB(I,ID4+ID) XPALB(I,ID9+ID) = XPALB(I,ID6+ID) 242 CONTINUE * end do ENDIF DO 244 ID = 1,IDIM FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFN*XPALB(I,ID3+ID) & + XPALB(I,ID6+ID) 244 CONTINUE * --- choc elementaire point_cercle_mobile * avec ou sans amortissement * on neglige la rotation (torsion) du cercle ELSE IF (ITYP.EQ.33 .OR. ITYP.EQ.34) THEN NPOA = IPLIB(I,1) NPOB = IPLIB(I,2) IGP = IPALB(I,2) IDIM = IPALB(I,3) INTER = IPALB(I,4) IF (ITYP.EQ.33) THEN ID1 = 6 ELSE ID1 = 7 ENDIF ID2 = ID1 + IDIM ID3 = ID1 + 2*IDIM ID4 = ID1 + 3*IDIM ID5 = ID1 + 4*IDIM ID6 = ID1 + 5*IDIM ID7 = ID1 + 6*IDIM ID8 = ID1 + 7*IDIM ID9 = ID1 + 8*IDIM * si pas d'adherence (c.a.d. glissement ou pas de contact) au pas precedent, * reactualisation de la position ecart origine d'adherence IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 330 ID=1,IDIM XPALB(I,ID7+ID) = (XPTB(NPOA,IND,ID) + & FEXPSM(NPOA,NPA,IND1,ID) ) & - ( XPTB(NPOB,IND,ID) + & FEXPSM(NPOB,NPA,IND1,ID)) 330 CONTINUE * end do ENDIF * calcul du deplacement relatif (PS) sur la normale au plan de section droite PS = 0.D0 DO 332 ID = 1,IDIM IDD2 = 3 + ID XDE2 = XPTB(NPOA,IND,ID) - XPTB(NPOB,IND,ID) XDE2 = XDE2 + FEXPSM(NPOa,NPA,IND1,ID) & - FEXPSM(NPOb,NPA,IND1,ID) XVALB(I,IND,IDD2) = XDE2 XDM2 = XPTB(NPOA,IND2,ID) - xptb (npob,ind2,id) XDM2 = XDM2 + FEXPSM(NPOA,NPAM1,INDM1,ID) & - FEXPSM(NPOb,NPAM1,INDM1,ID) XPALB(I,ID4+ID) = XDE2 - XDM2 XPALB(I,ID5+ID) = XDE2 - XPALB(I,ID7+ID) PS = PS + XDE2 * XPALB(I,ID1+ID) 332 CONTINUE * end do * calcul de la valeur du deplacement relatif suivant la normale * au plan tangent du contact (on tient compte de l'excentrement) XPME = 0.D0 PSXPME = 0.D0 DO 334 ID = 1,IDIM XDE2 = XPTB(NPOA,IND,ID) - XPTB(NPOB,IND,ID) XDE2 = XDE2 + FEXPSM(NPOa,NPA,IND1,ID) & - FEXPSM(NPOb,NPA,IND1,ID) XPRIM = XDE2 - PS * XPALB(I,ID1+ID) XPME = XPRIM - (XPALB(I,ID2+ID) - & XPALB(I,ID2+ID)*XPALB(I,ID1+ID)) XPALB(I,ID3+ID) = XPME PSXPME = PSXPME + XPME * XPME 334 CONTINUE * end do XDEP = SQRT(PSXPME) * normale au plan tangent * IF(XDEP.GT.1.D-20) THEN IF(XDEP.GT.xpetit) THEN DO 336 ID = 1,IDIM XPALB(I,ID3+ID) = XPALB(I,ID3+ID) / XDEP 336 CONTINUE * end do ENDIF * calcul de la vitesse relative tangentielle par derivee a gauche PSN = 0.D0 PSN0 = 0.D0 DO 338 ID = 1,IDIM PSN = PSN + XPALB(I,ID4+ID) * XPALB(I,ID3+ID) PSN0 = PSN0 + XPALB(I,ID5+ID) * XPALB(I,ID3+ID) 338 CONTINUE * end do DO 340 ID = 1,IDIM XPALB(I,ID4+ID) = (XPALB(I,ID4+ID) - & PSN * XPALB(I,ID3+ID))/PDTS2 XPALB(I,ID5+ID) = XPALB(I,ID5+ID) - & PSN0 * XPALB(I,ID3+ID) 340 CONTINUE * end do * calcul de la force de choc IF (ITYP.EQ.33) THEN CALL DYCHE3(XDEP,IDIM,IGP,XPALB,NLIAB,I,INTER & ,XFN,XFT,XPUS,iannul) ELSE XVITN = PSN / PDTS2 XVALB(I,IND,3) = XVITN CALL DYCHA3(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER & ,XFN,XFT,XPUS,iannul) ENDIF XVALB(I,IND,1) = XFN XVALB(I,IND,10) = ABS(XFT) XVALB(I,IND,12) = XPUS IPALB(I,2) = IGP IF (IGP .EQ. 1) THEN PS = 0.D0 DO 30 ID = 1,IDIM PS = PS + (XPALB(I,ID4+ID)*XPALB(I,ID4+ID)) 30 CONTINUE XVITT = SQRT(PS) ELSE XVITT = 0.D0 ENDIF XVALB(I,IND,11) = XVITT * si pas d'adherence (c.a.d. glissement ou pas de contact) * memorisation de la vitesse tangentielle et de la force tangentielle IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 342 ID = 1,IDIM XPALB(I,ID8+ID) = XPALB(I,ID4+ID) XPALB(I,ID9+ID) = XPALB(I,ID6+ID) 342 CONTINUE * end do ENDIF DO 344 ID = 1,IDIM FTOTB(NPOA,ID) = FTOTB(NPOA,ID) + XFN*XPALB(I,ID3+ID) & + XPALB(I,ID6+ID) FTOTB(NPOB,ID) = FTOTB(NPOB,ID) - XFN*XPALB(I,ID3+ID) & - XPALB(I,ID6+ID) 344 CONTINUE * * --- choc elementaire cercle_cercle_frottement * avec ou sans amortissement * ELSE IF (ITYP.EQ.25 .OR. ITYP.EQ.26) THEN NPOI = IPLIB(I,1) IGP = IPALB(I,2) IDIM = IPALB(I,3) INTER = IPALB(I,4) IF (ITYP.EQ.25) THEN ID1 = 6 ELSE ID1 = 7 ENDIF ID2 = ID1 + IDIM ID3 = ID1 + 2*IDIM ID4 = ID1 + 3*IDIM ID5 = ID1 + 4*IDIM ID6 = ID1 + 5*IDIM ID7 = ID1 + 6*IDIM ID8 = ID1 + 7*IDIM ID9 = ID1 + 8*IDIM ID10 = ID1 + 9*IDIM XRAYT = XPALB(I,ID10+1) XREXT = XPALB(I,2) * calcul du deplacement du point fibre neutre dans le plan du cercle * recuperation de la normale de choc au pas precedent DO 249 ID = 1,IDIM XN2(ID)= XPALB(I,ID3+ID) 249 CONTINUE * calcul de la normale de choc au pas courant PSXPN = 0.D0 DO 250 ID = 1,IDIM PSXPN = PSXPN + ( XPTB(NPOI,IND,ID) * XPALB(I,ID1+ID) ) 250 CONTINUE * end do PSXPME = 0.D0 DO 254 ID = 1,IDIM XXPME = ( XPTB(NPOI,IND,ID) - ( PSXPN * XPALB(I,ID1+ID) ) ) & - XPALB(I,ID2+ID) XPALB(I,ID3+ID) = XXPME PSXPME = PSXPME + ( XXPME * XXPME ) 254 CONTINUE * end do PSXPME = SQRT(PSXPME) * IF (PSXPME.GT.1D-20) THEN IF (PSXPME.GT.xpetit) THEN DO 256 ID = 1,IDIM XPALB(I,ID3+ID) = XPALB(I,ID3+ID) / PSXPME 256 CONTINUE * end do ENDIF * valeur du deplacement du point de contact suivant la normale XDEP = PSXPME + XRAYT * calcul du deplacement du point de contact au pas courant XPTP2(1) = XPTB(NPOI,IND,1) + (XPALB(I,ID3+1)*XRAYT) XPTP2(2) = XPTB(NPOI,IND,2) + (XPALB(I,ID3+2)*XRAYT) XPTP2(3) = XPTB(NPOI,IND,3) + (XPALB(I,ID3+3)*XRAYT) * calcul du deplacement du point de contact au pas precedent XPTPM2(1) = XPTB(NPOI,IND2,1) + (XN2(1)*XRAYT) XPTPM2(2) = XPTB(NPOI,IND2,2) + (XN2(2)*XRAYT) XPTPM2(3) = XPTB(NPOI,IND2,3) + (XN2(3)*XRAYT) * Vitesse *(-1) du point de contact appartenant a la structure mobile * due a la rotation absolue XVPC0(1) = (1.D0/ PDTS2) * & ( ( XPTB(NPOI,IND2,5) * XPALB(I,ID3+3) * XRAYT ) - & ( XPTB(NPOI,IND2,6) * XPALB(I,ID3+2) * XRAYT ) - & ( XPTB(NPOI,IND ,5) * XPALB(I,ID3+3) * XRAYT ) + & ( XPTB(NPOI,IND ,6) * XPALB(I,ID3+2) * XRAYT ) ) XVPC0(2) = (1.D0/ PDTS2) * & ( ( XPTB(NPOI,IND2,6) * XPALB(I,ID3+1) * XRAYT ) - & ( XPTB(NPOI,IND2,4) * XPALB(I,ID3+3) * XRAYT ) - & ( XPTB(NPOI,IND ,6) * XPALB(I,ID3+1) * XRAYT ) + & ( XPTB(NPOI,IND ,4) * XPALB(I,ID3+3) * XRAYT ) ) XVPC0(3) = (1.D0/ PDTS2) * & ( ( XPTB(NPOI,IND2,4) * XPALB(I,ID3+2) * XRAYT ) - & ( XPTB(NPOI,IND2,5) * XPALB(I,ID3+1) * XRAYT ) - & ( XPTB(NPOI,IND ,4) * XPALB(I,ID3+2) * XRAYT ) + & ( XPTB(NPOI,IND ,5) * XPALB(I,ID3+1) * XRAYT ) ) * si glissement au pas precedent, reactualisation de la position * origine d'adherence a l'aide du point de contact IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 257 ID=1,IDIM XPALB(I,ID7+ID) = XPTP2(ID) 257 CONTINUE ELSE PXVPC0 = 0.D0 DO 258 ID = 1,IDIM PXVPC0 = PXVPC0 + XVPC0(ID)*XPALB(I,ID1+ID) 258 CONTINUE DO 259 ID = 1,IDIM XVPCT(ID) = XVPC0(ID) - PXVPC0 * XPALB(I,ID1+ID) 259 CONTINUE XPALB(I,ID7+1) = XPALB(I,ID7+1) + & (XVPCT(1)*PDTS2)*(XREXT/(XREXT-XRAYT)) XPALB(I,ID7+2) = XPALB(I,ID7+2) + & (XVPCT(2)*PDTS2)*(XREXT/(XREXT-XRAYT)) XPALB(I,ID7+3) = XPALB(I,ID7+3) + & (XVPCT(3)*PDTS2)*(XREXT/(XREXT-XRAYT)) ENDIF * calcul du deplacement sur la normale au plan de section droite * et de l'ecart a la position orgine adherencee DO 260 ID = 1,IDIM IDD1 = 3 + ID IDD2 = 6 + ID IDD3 = 15 + ID XVALB(I,IND,IDD1) = XPTB(NPOI,IND,ID) XVALB(I,IND,IDD2) = (XPTP2(ID) - XPTPM2(ID) ) / PDTS2 XVALB(I,IND,IDD3) = XPTB(NPOI,IND,ID+3) XPALB(I,ID4+ID) = XPTP2(ID) - XPTPM2(ID) & - (XVPC0(ID) * PDTS2) XPALB(I,ID5+ID) = XPTP2(ID) - XPALB(I,ID7+ID) 260 CONTINUE * end do * calcul de la vitesse tangentielle par derivee a gauche * et de l'ecart a la position orgine adherencee PSN = 0.D0 PSN0 = 0.D0 DO 262 ID = 1,IDIM PSN = PSN + XPALB(I,ID4+ID) * XPALB(I,ID3+ID) PSN0 = PSN0 + XPALB(I,ID5+ID) * XPALB(I,ID3+ID) 262 CONTINUE * end do DO 264 ID = 1,IDIM XPALB(I,ID4+ID) = (XPALB(I,ID4+ID) - & PSN * XPALB(I,ID3+ID))/PDTS2 XPALB(I,ID5+ID) = XPALB(I,ID5+ID) - & PSN0 * XPALB(I,ID3+ID) 264 CONTINUE * end do * calcul de la force de choc IF (ITYP.EQ.25) THEN CALL DYCHE3(XDEP,IDIM,IGP,XPALB,NLIAB,I,INTER & ,XFN,XFT,XPUS,iannul) ELSE XVITN = PSN / PDTS2 XVALB(I,IND,3) = XVITN CALL DYCHA3(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,INTER & ,XFN,XFT,XPUS,iannul) ENDIF XVALB(I,IND,1) = XFN XVALB(I,IND,10) = ABS(XFT) XVALB(I,IND,12) = XPUS IPALB(I,2) = IGP * si glissement, memorisation de la vitesse tangentielle et de la force * tangentielle IF (IGP.EQ.1) THEN DO 266 ID = 1,IDIM XPALB(I,ID8+ID) = XPALB(I,ID4+ID) XPALB(I,ID9+ID) = XPALB(I,ID6+ID) 266 CONTINUE ELSE DO 267 ID = 1,IDIM XPALB(I,ID9+ID) = 0.D0 267 CONTINUE * end do ENDIF DO 268 ID = 1,IDIM XFOR = ( XFN * XPALB(I,ID3+ID) ) + XPALB(I,ID6+ID) FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFOR XFNT(ID) = XPALB (I ,ID6+ID) 268 CONTINUE * end do XAPP1 = XRAYT * XPALB(I,ID3+1) XAPP2 = XRAYT * XPALB(I,ID3+2) XAPP3 = XRAYT * XPALB(I,ID3+3) XAPFP1 = ( XAPP2 * XFNT(3) ) - ( XAPP3 * XFNT(2) ) XAPFP2 = ( XAPP3 * XFNT(1) ) - ( XAPP1 * XFNT(3) ) XAPFP3 = ( XAPP1 * XFNT(2) ) - ( XAPP2 * XFNT(1) ) XVALB(I,IND,13) = XAPFP1 XVALB(I,IND,14) = XAPFP2 XVALB(I,IND,15) = XAPFP3 FTOTB(NPOI,4) = FTOTB(NPOI,4) + XAPFP1 FTOTB(NPOI,5) = FTOTB(NPOI,5) + XAPFP2 FTOTB(NPOI,6) = FTOTB(NPOI,6) + XAPFP3 * * --- choc ........... * * else if (ityp.eq. ) then * ....... * ....... * ENDIF * END