C DEVFB2 SOURCE PV 08/04/18 21:15:01 6094 SUBROUTINE DEVFB2(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB, & NPLB,IND,IND1,INDM1,NPA,NPAM1,IND2,PDT,PDTS2, & FEXPSM,NPC1,XABSCI,XORDON,NIP,I,iannul) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) *--------------------------------------------------------------------* * * * Op{rateur DYNE : algorithme de Fu - de Vogelaere * * ________________________________________________ * * * * Calcul des forces de choc sur base B pour les liaisons de * * type POINT_POINT. * * * * Param}tres: * * * * e ITYP type de la liaison. * * es FTOTB Forces ext{rieures totalis{es sur la base B. * * e XPTB Tableau des d{placements des points * * e IPALB Renseigne sur la liaison. * * e IPLIB Tableau contenant les num{ros "DYNE" de la liaison. * * e XPALB Tableau contenant les param}tres 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 cr{ation: * * * * Lionel VIVAN : le 22 Septembre 1989 : Cr{ation * * * *--------------------------------------------------------------------* * 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 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP) * * * --- choc {l{mentaire POINT_POINT avec amortissement * IF (ITYP.EQ.11 .OR. ITYP.EQ.111) THEN NPOA = IPLIB(I,1) NPOB = IPLIB(I,2) IDIM = IPALB(I,3) IPERM = IPALB(I,4) XRAID = XPALB(I,1) XJEU = XPALB(I,2) XAMO = XPALB(I,3) XDEP = 0.D0 XDEPM1 = 0.D0 DO 20 ID = 1,IDIM IDA = 3 + ID IDB = 3 + IDIM + ID XDEA = XPTB(NPOA,IND,ID) XDMA = XPTB(NPOA,IND2,ID) XDEB = XPTB(NPOB,IND,ID) XDMB = XPTB(NPOB,IND2,ID) XVALB(I,IND,IDA) = XDEA XVALB(I,IND,IDB) = XDEB XDEA = XDEA + FEXPSM(NPOA,NPA,IND1,ID) XDEB = XDEB + FEXPSM(NPOB,NPA,IND1,ID) XDMA = XDMA + FEXPSM(NPOA,NPAM1,INDM1,ID) XDMB = XDMB + FEXPSM(NPOB,NPAM1,INDM1,ID) XDEP = XDEP + (XDEA - XDEB) * XPALB(I,3+ID) XDEPM1 = XDEPM1 + (XDMA - XDMB) * XPALB(I,3+ID) 20 CONTINUE * end do XVIT = (XDEP - XDEPM1) / PDTS2 XVALB(I,IND,3) = XVIT IF (ITYP.EQ.11) THEN CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFLA,IPERM,iannul) ELSE CALL DYCHAM2(XDEP,XVIT,XJEU,XAMO,XABSCI,XORDON,NIP, & NLIAB,I, XFLA,IPERM,iannul) ENDIF XFLB = -1.D0 * XFLA XVALB(I,IND,1) = XFLA XVALB(I,IND,2) = XFLB DO 22 ID = 1,IDIM FTOTB(NPOA,ID) = FTOTB(NPOA,ID) + XFLA * XPALB(I,3+ID) FTOTB(NPOB,ID) = FTOTB(NPOB,ID) + XFLB * XPALB(I,3+ID) 22 CONTINUE * end do * * --- choc {l{mentaire POINT_POINT_FROTTEMENT avec ou sans amortissement * ELSEIF ((ITYP.EQ.13) .or. ( ityp .eq. 113)) THEN NPOA = IPLIB(I,1) NPOB = IPLIB(I,2) IGP = IPALB(I,2) IDIM = IPALB(I,3) ID1 = 7 ID2 = ID1 + IDIM ID3 = ID1 + 2*IDIM ID4 = ID1 + 3*IDIM ID5 = ID1 + 4*IDIM ID6 = ID1 + 5*IDIM ID7 = ID1 + 6*IDIM * Si glissement au pas pr{c{dent, r{actualisation de la position-ecart * origine d'adh{rence IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 30 ID=1,IDIM XPALB(I,ID5+ID) =(XPTB(NPOa,IND,ID) + & FEXPSM(NPOa,NPA,IND1,ID) ) & - ( XPTB(NPOb,IND,ID) + & FEXPSM(NPOb,NPA,IND1,ID)) 30 CONTINUE * end do ENDIF * Calcul de l'enfoncement relatif et de la vitesse normale relative XDEP = 0.D0 PSN = 0.D0 PSN0 = 0.D0 DO 32 ID = 1,IDIM IDD1 = 3 + ID xvalb(i,ind,idd1) = XPTB(NPOa,IND,ID) xvalb(i,ind,idd1 + idim) = XPTB(NPOb,IND,ID) XDE2 = XPTB(NPOa,IND,ID) - xptb (npob,ind,id) XDE2 = XDE2 + FEXPSM(NPOa,NPA,IND1,ID) & - FEXPSM(NPOb,NPA,IND1,ID) XDm2 = XPTB(NPOa,IND2,ID) - xptb (npob,ind2,id) XDM2 = XDM2 + FEXPSM(NPOA,NPAM1,INDM1,ID) & - FEXPSM(NPOb,NPAM1,INDM1,ID) XPALB(I,ID2+ID) = XDE2 - XDM2 XPALB(I,ID3+ID) = XDE2 - XPALB(I,ID5+ID) XDEP = XDEP + XDE2 * XPALB(I,ID1+ID) PSN = PSN + XPALB(I,ID2+ID) * XPALB(I,ID1+ID) PSN0 = PSN0 + XPALB(I,ID3+ID) * XPALB(I,ID1+ID) 32 CONTINUE * end do * Projette la vitesse relative * et la variation de d{placement relatif par rapport a * l' ecart origine d'adh{rence sur le plan tangent DO 34 ID = 1,IDIM XPALB(I,ID2+ID) = (XPALB(I,ID2+ID) - PSN * XPALB(I,ID1+ID)) & / PDTS2 XPALB(I,ID3+ID) = XPALB(I,ID3+ID) - PSN0 * XPALB(I,ID1+ID) 34 CONTINUE * end do XVITN = PSN / PDTS2 XVALB(I,IND,3) = XVITN IF (ITYP.EQ.13) THEN CALL DYCHA4(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS & ,iannul) ELSE CALL DYCHA41(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS, & XABSCI,XORDON,NIP,iannul) ENDIF Xfla = XFN Xflb = -1d0 * XFN XVALB(I,IND,1) = XFla XVALB(I,IND,2) = XFlb XVALB(I,IND,10) = ABS(XFT) XVALB(I,IND,12) = XPUS IPALB(I,2) = IGP * Si glissement, m{morisation de la vitesse tangentielle et de la force * tangentielle IF (IGP.EQ.1.OR.IGP.EQ.-1) THEN DO 36 ID = 1,IDIM XPALB(I,ID6+ID) = XPALB(I,ID2+ID) XPALB(I,ID7+ID) = XPALB(I,ID4+ID) 36 CONTINUE * end do ENDIF DO 38 ID = 1,IDIM FTOTB(NPOa,ID) = FTOTB(NPOa,ID) + XFla* XPALB(I,ID1+ID) & + XPALB(I,ID4+ID) FTOTB(NPOb,ID) = FTOTB(NPOb,ID) + XFlb* XPALB(I,ID1+ID) & - XPALB(I,ID4+ID) 38 CONTINUE ***************************************************************** * ** modele de Nedjai. On ne garantit rien ELSEIF (ITYP.EQ.-13) THEN NPOA = IPLIB(I,1) NPOB = IPLIB(I,2) IGP = IPALB(I,2) IDIM = IPALB(I,3) ID1 = 7 ID2 = ID1 + IDIM ID3 = ID1 + 2*IDIM ID4 = ID1 + 3*IDIM ID5 = ID1 + 4*IDIM ID6 = ID1 + 5*IDIM ID7 = ID1 + 6*IDIM * Si glissement au pas pr{c{dent, r{actualisation de la position-ecart * origine d'adh{rence IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 130 ID=1,IDIM XPALB(I,ID5+ID) =(XPTB(NPOa,IND,ID) + & FEXPSM(NPOa,NPA,IND1,ID) ) & - ( XPTB(NPOb,IND,ID) + & FEXPSM(NPOb,NPA,IND1,ID)) 130 CONTINUE * end do ENDIF * Calcul de l'enfoncement relatif et de la vitesse normale relative XDEP = 0.D0 PSN = 0.D0 PSN0 = 0.D0 DO 132 ID = 1,IDIM IDD1 = 3 + ID xvalb(i,ind,idd1) = XPTB(NPOa,IND,ID) xvalb(i,ind,idd1 + idim) = XPTB(NPOb,IND,ID) XDE2 = XPTB(NPOa,IND,ID) - xptb (npob,ind,id) XDE2 = XDE2 + FEXPSM(NPOa,NPA,IND1,ID) & - FEXPSM(NPOb,NPA,IND1,ID) XDm2 = XPTB(NPOa,IND2,ID) - xptb (npob,ind2,id) XDM2 = XDM2 + FEXPSM(NPOA,NPAM1,INDM1,ID) & - FEXPSM(NPOb,NPAM1,INDM1,ID) *** * WRITE(6,*) 'devfb2 XDE2 ', XDE2 * WRITE(6,*) 'devfb2 XDM2 ', XDM2 *** XPALB(I,ID2+ID) = XDE2 - XDM2 *A ENLEVER**XPALB(I,ID3+ID) = XDE2 - XPALB(I,ID5+ID) XDEP = XDEP + XDE2 * XPALB(I,ID1+ID) PSN = PSN + XPALB(I,ID2+ID) * XPALB(I,ID1+ID) PSN0 = PSN0 + XPALB(I,ID3+ID) * XPALB(I,ID1+ID) 132 CONTINUE * end do * Projette la vitesse relative * et la variation de d{placement relatif par rapport a * l' ecart origine d'adh{rence sur le plan tangent DO 134 ID = 1,IDIM XPALB(I,ID2+ID) = (XPALB(I,ID2+ID) - PSN * XPALB(I,ID1+ID)) & / PDTS2 *A enlever**XPALB(I,ID3+ID) = XPALB(I,ID3+ID) - PSN0 * XPALB(I,ID1+ID) **** RAJOUT D UN ORDRE D IMPRESSION ************************************* * WRITE(6,*) 'devfb2 XPALB',XPALB(I,ID1+ID) * WRITE(6,*) 'devfb2 XPALB',XPALB(I,ID2+ID) * WRITE(6,*) 'devfb2 XPALB',XPALB(I,ID3+ID) 134 CONTINUE * end do XVITN = PSN / PDTS2 XVALB(I,IND,3) = XVITN ***RAJOUT DE XPAS DANS CALL DYCHA4 CALL DGCHA4(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS & ,iannul,PDTS2) Xfla = XFN Xflb = -1d0 * XFN XVALB(I,IND,1) = XFla XVALB(I,IND,2) = XFlb XVALB(I,IND,10) = ABS(XFT) XVALB(I,IND,12) = XPUS IPALB(I,2) = IGP * Si glissement, m{morisation de la vitesse tangentielle et de la force * tangentielle *AENLEV**IF (IGP.EQ.1.OR.IGP.EQ.-1) THEN DO 136 ID = 1,IDIM **ON A REMPLACE ID6 PAR ID3 (NEDJAI) XPALB(I,ID3+ID) = XPALB(I,ID2+ID) XPALB(I,ID7+ID) = XPALB(I,ID4+ID) 136 CONTINUE * end do **AENLEV**ENDIF DO 138 ID = 1,IDIM FTOTB(NPOa,ID) = FTOTB(NPOa,ID) + XFla* XPALB(I,ID1+ID) & + XPALB(I,ID4+ID) FTOTB(NPOb,ID) = FTOTB(NPOb,ID) + XFlb* XPALB(I,ID1+ID) & - XPALB(I,ID4+ID) 138 CONTINUE * * --- choc ........... * * ELSE IF (ITYP.EQ. ) THEN * ....... * ....... * ENDIF * END