C DEVFB1 SOURCE CHAT 05/01/12 22:45:16 5004 SUBROUTINE DEVFB1(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_PLAN. * * * * 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 * * Bertrand BEAUFILS : le 31 Mai 1990 : Ajout frottement sec* * * *--------------------------------------------------------------------* * 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),XPTPP2(3),XFNT(3),XPSMP(3),XPSMPM(3) REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP) * * --- choc {l{mentaire POINT_PLAN * avec ou sans amortissement * IF (ITYP.EQ.1 .OR. ITYP.EQ.100 & .OR. ITYP.EQ.101 .OR. ITYP.EQ.102) THEN NPOI = IPLIB(I,1) IDIM = IPALB(I,3) IPERM = IPALB(I,4) XRAID = XPALB(I,1) XJEU = XPALB(I,2) XAMO = XPALB(I,3) IF (ITYP.EQ.1 .OR. ITYP.EQ.102) THEN ID1 = 3 ** ianis ELSE ID1 = 4 XSEUIL = XPALB(I,ID1) XDPLAS = XPALB(I,(ID1+IDIM+1)) ENDIF * XDEP = 0.D0 DO 10 ID = 1,IDIM IDD1 = 3 + ID XDE2 = XPTB(NPOI,IND,ID) XVALB(I,IND,IDD1) = XDE2 XDE2 = XDE2 + FEXPSM(NPOI,NPA,IND1,ID) XDEP = XDEP + XDE2 * XPALB(I,ID1+ID) 10 CONTINUE * end do XDEPM1 = 0.D0 DO 12 ID = 1,IDIM XDM2 = XPTB(NPOI,IND2,ID) + FEXPSM(NPOI,NPAM1,INDM1,ID) XDEPM1 = XDEPM1 + XDM2 * XPALB(I,ID1+ID) 12 CONTINUE * end do XVIT = (XDEP - XDEPM1) / PDTS2 XVALB(I,IND,3) = XVIT ** ianis IF (ITYP.EQ.1) THEN CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul) ELSE IF (ITYP.EQ.102) THEN CALL DYCHAM2(XDEP,XVIT,XJEU,XAMO,XABSCI,XORDON,NIP, & NLIAB,I,XFL,IPERM,iannul) ELSE IF (ITYP.EQ.101) THEN CALL DYCHPL2(XDEP,XVIT,XDPLAS,XRAID,XJEU,XAMO,XSEUIL, & XABSCI,XORDON,NIP,NLIAB,I,XFL,IPERM,iannul) ELSE IF (ITYP.EQ.100) THEN CALL DYCHPL(XDEP,XVIT,XDPLAS,XRAID,XJEU,XAMO,XSEUIL, & XFL,IPERM,iannul) ENDIF IF (ITYP.EQ.100 .OR. ITYP.EQ.101) THEN XPALB(I,(ID1+IDIM+1)) = XDPLAS XVALB(I,IND,13) = XDPLAS ENDIF * XVALB(I,IND,1) = XFL DO 14 ID = 1,IDIM FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFL * XPALB(I,ID1+ID) 14 CONTINUE * end do * * --- choc {l{mentaire POINT_PLAN_FROTTEMENT * avec ou sans amortissement * ELSE IF (ITYP.EQ.3 .OR. ITYP.EQ.103) THEN NPOI = IPLIB(I,1) 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 * origine d'adh{rence IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 30 ID=1,IDIM XPALB(I,ID5+ID) = XPTB(NPOI,IND,ID) + & FEXPSM(NPOI,NPA,IND1,ID) 30 CONTINUE * end do ENDIF * Calcul de l'enfoncement et de la vitesse normale XDEP = 0.D0 PSN = 0.D0 PSN0 = 0.D0 DO 32 ID = 1,IDIM IDD1 = 3 + ID XDE2 = XPTB(NPOI,IND,ID) XDM2 = XPTB(NPOI,IND2,ID) XVALB(I,IND,IDD1) = XDE2 XDE2 = XDE2 + FEXPSM(NPOI,NPA,IND1,ID) XDM2 = XDM2 + FEXPSM(NPOI,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 et la variation de d{placement par rapport a * la position 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.3) 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 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,ID2+ID)*XPALB(I,ID2+ID)) 20 CONTINUE XVITT = SQRT(PS) ELSE XVITT = 0.D0 ENDIF XVALB(I,IND,11) = XVITT * 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(NPOI,ID) = FTOTB(NPOI,ID) + XFN * XPALB(I,ID1+ID) & + XPALB(I,ID4+ID) 38 CONTINUE * end do * * --- choc {l{mentaire CERCLE_PLAN_FROTTEMENT * avec ou sans amortissement * ELSE IF (ITYP.EQ.5 .OR. ITYP.EQ.6) THEN NPOI = IPLIB(I,1) IGP = IPALB(I,2) IDIM = IPALB(I,3) IF (ITYP.EQ.5) 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 + 1 + 7*IDIM * calcul du d{placement du point de contact au pas pr{c{dent XMP1 = XPALB(I,ID1+1) * XPALB(I,ID8) XMP2 = XPALB(I,ID1+2) * XPALB(I,ID8) XMP3 = XPALB(I,ID1+3) * XPALB(I,ID8) XPTPP2(1) = XPTB(NPOI,IND2,1) + & ( ( XPTB(NPOI,IND2,5) * XMP3 ) - & ( XPTB(NPOI,IND2,6) * XMP2 ) ) XPTPP2(2) = XPTB(NPOI,IND2,2) + & ( ( XPTB(NPOI,IND2,6) * XMP1 ) - & ( XPTB(NPOI,IND2,4) * XMP3 ) ) XPTPP2(3) = XPTB(NPOI,IND2,3) + & ( ( XPTB(NPOI,IND2,4) * XMP2 ) - & ( XPTB(NPOI,IND2,5) * XMP1 ) ) XPSMPM(1) = FEXPSM(NPOI,NPAM1,INDM1,1) + & ( ( FEXPSM(NPOI,NPAM1,INDM1,5) * XMP3 ) - & ( FEXPSM(NPOI,NPAM1,INDM1,6) * XMP2 ) ) XPSMPM(2) = FEXPSM(NPOI,NPAM1,INDM1,2) + & ( ( FEXPSM(NPOI,NPAM1,INDM1,6) * XMP1 ) - & ( FEXPSM(NPOI,NPAM1,INDM1,4) * XMP3 ) ) XPSMPM(3) = FEXPSM(NPOI,NPAM1,INDM1,3) + & ( ( FEXPSM(NPOI,NPAM1,INDM1,4) * XMP2 ) - & ( FEXPSM(NPOI,NPAM1,INDM1,5) * XMP1 ) ) * calcul du d{placement du point de contact au pas courant XPTP2(1) = XPTB(NPOI,IND,1) + & ( ( XPTB(NPOI,IND,5) * XMP3 ) - & ( XPTB(NPOI,IND,6) * XMP2 ) ) XPTP2(2) = XPTB(NPOI,IND,2) + & ( ( XPTB(NPOI,IND,6) * XMP1 ) - & ( XPTB(NPOI,IND,4) * XMP3 ) ) XPTP2(3) = XPTB(NPOI,IND,3) + & ( ( XPTB(NPOI,IND,4) * XMP2 ) - & ( XPTB(NPOI,IND,5) * XMP1 ) ) XPSMP(1) = FEXPSM(NPOI,NPA,IND1,1) + & ( ( FEXPSM(NPOI,NPA,IND1,5) * XMP3 ) - & ( FEXPSM(NPOI,NPA,IND1,6) * XMP2 ) ) XPSMP(2) = FEXPSM(NPOI,NPA,IND1,2) + & ( ( FEXPSM(NPOI,NPA,IND1,6) * XMP1 ) - & ( FEXPSM(NPOI,NPA,IND1,4) * XMP3 ) ) XPSMP(3) = FEXPSM(NPOI,NPA,IND1,3) + & ( ( FEXPSM(NPOI,NPA,IND1,4) * XMP2 ) - & ( FEXPSM(NPOI,NPA,IND1,5) * XMP1 ) ) * * Si glissement au pas pr{c{dent, r{actualisation de la position * origine d'adh{rence IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN DO 50 ID=1,IDIM XPALB(I,ID5+ID) = XPTP2(ID) + XPSMP(ID) 50 CONTINUE * end do ENDIF * Calcul de l'enfoncement et de la vitesse normale XDEP = 0.D0 PSN = 0.D0 PSN0 = 0.D0 DO 52 ID = 1,IDIM IDD1 = 3 + ID IDD2 = 6 + ID IDD3 = 15 + ID XDE2 = XPTP2(ID) XDM2 = XPTPP2(ID) XVALB(I,IND,IDD1) = XPTB(NPOI,IND,ID) XVALB(I,IND,IDD2) = (xDE2 - xdm2 ) / pdts2 XVALB(I,IND,IDD3) = XPTB(NPOI,IND,ID+3) XDE2 = XDE2 + XPSMP(ID) XDM2 = XDM2 + XPSMPM(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) 52 CONTINUE * end do * Projette la vitesse et la variation de d{placement par rapport a * la position origine d'adh{rence sur le plan tangent DO 54 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) 54 CONTINUE * end do IF (ITYP.EQ.5) THEN CALL DYCHE4(XDEP,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS, & iannul) ELSE XVITN = PSN / PDTS2 XVALB(I,IND,3) = XVITN CALL DYCHA4(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,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, m{morisation de la vitesse tangentielle et de la force * tangentielle IF (IGP.EQ.1.OR.IGP.EQ.-1) THEN DO 56 ID = 1,IDIM XPALB(I,ID6+ID) = XPALB(I,ID2+ID) XPALB(I,ID7+ID) = XPALB(I,ID4+ID) 56 CONTINUE * end do ENDIF DO 58 ID = 1,IDIM XFOR = ( XFN * XPALB(I,ID1+ID) ) + XPALB(I,ID4+ID) XFNT(ID) = XFOR FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFOR 58 CONTINUE * end do XMPFP1 = ( XMP2 * XFNT(3) ) - ( XMP3 * XFNT(2) ) XMPFP2 = ( XMP3 * XFNT(1) ) - ( XMP1 * XFNT(3) ) XMPFP3 = ( XMP1 * XFNT(2) ) - ( XMP2 * XFNT(1) ) XVALB(I,IND,13) = XMPFP1 XVALB(I,IND,14) = XMPFP2 XVALB(I,IND,15) = XMPFP3 FTOTB(NPOI,4) = FTOTB(NPOI,4) + XMPFP1 FTOTB(NPOI,5) = FTOTB(NPOI,5) + XMPFP2 FTOTB(NPOI,6) = FTOTB(NPOI,6) + XMPFP3 * * --- choc {l{mentaire POINT_PLAN_FLUIDE * ELSE IF (ITYP.EQ.7) THEN NPOI = IPLIB(I,1) IDIM = IPALB(I,3) XINER = XPALB(I,1) XCONV = XPALB(I,2) XVISC = XPALB(I,3) XPCEL = XPALB(I,4) XPCRA = XPALB(I,5) XJEU = XPALB(I,6) ID1 = 6 ID2 = ID1 + IDIM XDEP = 0.D0 DO 70 ID = 1,IDIM XDE2 = XPTB(NPOI,IND,ID) XVALB(I,IND,ID) = XDE2 XDE2 = XDE2 + FEXPSM(NPOI,NPA,IND1,ID) XDEP = XDEP + XDE2 * XPALB(I,ID1+ID) 70 CONTINUE * end do XDEPM1 = XPALB(I,ID2+1) XVITM1 = XPALB(I,ID2+2) XACCM1 = XPALB(I,ID2+3) XVIT = (XDEP - XDEPM1) / PDTS2 XACC = (XVIT - XVITM1) / PDTS2 XPALB(I,ID2+1) = XDEP XPALB(I,ID2+2) = XVIT XPALB(I,ID2+3) = XACC XDH = XJEU - XDEP * Calcul de la force d'inertie CALL DYFINE(XDH,XDEP,XACC,XJEU,XINER,XFI,iannul) * Calcul de la force de convection CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFC,iannul) * Calcul de la force de viscosit{ CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFV,iannul) * Calcul de la force de perte de charge CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFP,iannul) * XFL = XFI + XFC + XFV + XFP DO 72 ID = 1,IDIM FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFL * XPALB(I,ID1+ID) 72 CONTINUE * end do XVALB(I,IND,IDIM+1) = XVIT XVALB(I,IND,IDIM+2) = XACC XVALB(I,IND,IDIM+3) = XFI XVALB(I,IND,IDIM+4) = XFC XVALB(I,IND,IDIM+5) = XFV XVALB(I,IND,IDIM+6) = XFP * * --- choc ........... * * ELSE IF (ITYP.EQ. ) THEN * ....... * ....... * ENDIF * END