C D2VLFA SOURCE CHAT 05/01/12 22:35:41 5004 C DEVLFA SOURCE IANIS 96/08/08 21:15:06 2249 SUBROUTINE D2VLFA(Q1,Q2,FTOTA,XM,NA1,IPALA,IPLIA,XPALA,XVALA, & NLIAA,XDT,NPAS,IND,FINERT,IVINIT,FTEST,FTOTA0) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) *--------------------------------------------------------------------* * * * Operateur DYNE : algorithme de Fu - de Vogelaere * * ________________________________________________ * * * * Calcul des forces de choc base A. * * * * Param}tres: * * * * e Q1(.,.) Vecteur des deplacements generalises. * * e Q2(.,.) Vecteur des vitesses generalisees. * * es FTOTA Forces ext{rieures totalis{es sur la base A. * * e XM Vecteur des masses g{n{ralis{es. * * e NA1 Nombre total d'inconnues en base A. * * e IPALA Renseigne sur le type de la liaison. * * e IPLIA Tableau contenant les num{ros "DYNE" de la liaison. * * e XPALA Tableau contenant les param}tres de la liaison. * * es XVALA Tableau contenant les variables internes des liaisons * * e NLIAA Nombre de liaisons sur la base A. * * e XDT Valeurs des pas de temps * * e NPAS Numero du pas de temps * * e IND Indice du demi-pas de temps * * = 2 si 1er demi-pas = 1 si 2eme demi-pas * * es FINERT Forces d'inertie * * e IVINIT =1 si vitesses initiales, =0 sinon * * * * Auteur, date de cr{ation: * * * * Lionel VIVAN, le 20 aout 1989. * * E. de LANGRE, 15 fevrier 91 correction jeu negatif * * Lionel VIVAN, 12 mars 1991 ajout de la liaison FLUIDE * * Denis ROBERT, 30 avril 1992 ajout liaison POLYNOMIALE * * * * remarque: * * ========= * * Si jeu n{gatif (cas particulier de la base A ou il n'y a pas de * * normale), on renverse les variables avec XNORM. * * * *--------------------------------------------------------------------* * INTEGER IPALA(NLIAA,*),IPLIA(NLIAA,*) REAL*8 XPALA(NLIAA,*),Q1(NA1,*),Q2(NA1,*),FTOTA(NA1,*),XM(*) REAL*8 XVALA(NLIAA,4,*),XDT(*),FINERT(NA1,*) PARAMETER (XZERO = 0.D0) REAL*8 FTest(nA1,4) REAL*8 ftOTa0(NA1,4) * * Boucle sur le nombre de liaisons * XFIN = 0.D0 PDT = XDT(NPAS) c PDTS2 = 0.5D0 * PDT pdts2 = pdt DO 10 I = 1,NLIAA ITYP = IPALA(I,1) icond= IPALa(I,2) iannul= 0 * --- cas des liaisons conditionnelles if (icond .eq. 1 ) then DO 101 j = 4,20 jliai = ipala(i,j) jpliai = abs ( jliai) if ( jliai . EQ. 0 ) then goto 101 else jtyp = ipala(jpliai,1) do 102 ik = 1,nA1 do 103 jk = 1,4 ftest(ik,jk) = 0d0 ftota0 (ik,jk) = ftota(ik,jk) 103 continue 102 continue * CALCUL DES LIAISONS POUR TEST * * ****** choc {l{mentaire POINT_PLAN sans amortissement * IF (JTYP.EQ.1) THEN XRAID = XPALA(jpliai,1) XJEU = XPALA(jpliai,2) XNORM = 1. IF (XJEU.LT.0D0) THEN XNORM = -1. XJEU = -1*XJEU ENDIF INA1 = IPLIA(jpliai,1) XDEP = XNORM*Q1(INA1,IND) CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul) XVALA(jpliai,IND,1) = XNORM*XFL XVALA(jpliai,IND,4) = XNORM*XDEP FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL * * ------ choc {l{mentaire POINT_PLAN avec amortissement * ELSE IF (JTYP.EQ.2) THEN XRAID = XPALA(jpliai,1) XJEU = XPALA(jpliai,2) XAMO = XPALA(jpliai,3) XNORM = 1. IPERM = 0 IF (XJEU.LT.XZERO) THEN XNORM = -1. XJEU = -1*XJEU ENDIF INA1 = IPLIA(jpliai,1) XDEP = XNORM * Q1(INA1,IND) IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN XVIT = XNORM * Q2(INA1,IND) ELSE IND2 = IND + 1 XDEPM1 = XNORM * Q1(INA1,IND2) XVIT = (XDEP - XDEPM1) / PDTS2 ENDIF XVALA(jpliai,IND,3) = XNORM*XVIT CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul) *+* IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0 IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0 *+* XVALA(jpliai,IND,1) = XNORM*XFL XVALA(jpliai,IND,4) = XNORM*XDEP FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL * * ------ choc {l{mentaire POINT_PLAN_FLUIDE * ELSE IF (JTYP.EQ.3) THEN XINER = XPALA(jpliai,1) XCONV = XPALA(jpliai,2) XVISC = XPALA(jpliai,3) XPCEL = XPALA(jpliai,4) XPCRA = XPALA(jpliai,5) XJEU = XPALA(jpliai,6) INA1 = IPLIA(jpliai,1) XDEP = Q1(INA1,IND) IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN XVIT = Q2(INA1,IND) ELSE IND2 = IND + 1 XDEPM1 = Q1(INA1,IND2) XVIT = (XDEP - XDEPM1) / PDTS2 ENDIF IF (XJEU.GT.0D0) THEN XDH= XJEU - XDEP XNORM = 1. ELSE XDH= XDEP - XJEU XNORM = -1. ENDIF * Calcul de la masse ajout{e XXIN = -XINER / XDH FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN * Calcul de la force de convection CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul) * Calcul de la force de viscosit{ CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul) * Calcul de la force de perte de charge CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,iannul) XFL = (XFCO* XNORM ) + XFVI + XFPE XVALA(jpliai,IND,1) = XDEP XVALA(jpliai,IND,2) = XVIT XVALA(jpliai,IND,3) = XXIN XVALA(jpliai,IND,4) = XFCO*XNORM XVALA(jpliai,IND,5) = XFVI XVALA(jpliai,IND,6) = XFPE FTest(INA1,IND) =Ftest(iNA1,IND) + XFL * * ------ force {l{mentaire de COUPLAGE EN VITESSE * ELSE IF (JTYP.EQ.4) THEN INA1 = IPLIA(jpliai,1) INA2 = IPLIA(jpliai,2) XDEP = Q1(INA2,IND) IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN XVIT = Q2(INA2,IND) ELSE IND2 = IND + 1 XDEPM1 = Q1(INA2,IND2) XVIT = (XDEP - XDEPM1) / PDTS2 ENDIF XCPLGE = XPALA(jpliai,1) CALL DYCPLV(XVIT,XCPLGE,XFL,iannul) XVALA(jpliai,IND,3) = XVIT XVALA(jpliai,IND,1) = XFL XVALA(jpliai,IND,4) = XDEP FTest(INA1,IND) =Ftest(INA1,IND) + XFL * * ------ force {l{mentaire de COUPLAGE EN DEPLACEMENT * ELSE IF (JTYP.EQ.5) THEN INA1 = IPLIA(jpliai,1) INA2 = IPLIA(jpliai,2) XDEP = Q1(INA2,IND) XCPLGE = XPALA(jpliai,1) CALL DYCPLD(XDEP,XCPLGE,XFL,iannul) XVALA(jpliai,IND,1) = XFL XVALA(jpliai,IND,4) = XDEP FTest(INA1,IND) =Ftest(INA1,IND) + XFL * * ------ force elementaire de type POLYNOMIALE * ELSE IF (JTYP.EQ.6) THEN * WRITE(IOIMP,*)'DEVLFA : IPLIA(',jpliai,'1)=',IPLIA(jpliai,1) * nombre de modes "origine" NMOD = IPALA(I,2) CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,XDT, & NPAS,jpliai,NMOD,FTest,IVINIT,iannul) * * ------ choc {l{mentaire ... * * ELSE IF (JTYP.EQ. ) THEN * ....... * ....... * ENDIF ***** tests pour voir si on annulera xff = 0d0 do 104 ik = 1,na1 do 105 jk = 1,4 xff = xff + ( ftest(ik,jk) ** 2) 105 continue 104 continue xff = xff ** .5 if ( ((xff .le. 1e-20 ) .and. ( jliai .gt. 0) ) & .OR. ((xff .gt. 1e-20 ) .and. ( jliai .lt. 0) ) ) & then iannul = 1 endif endif 101 continue * --- fin du cas des liaisons conditionnelles endif * * ------ choc {l{mentaire POINT_PLAN sans amortissement * IF (ITYP.EQ.1) THEN XRAID = XPALA(I,1) XJEU = XPALA(I,2) XNORM = 1. IF (XJEU.LT.0D0) THEN XNORM = -1. XJEU = -1*XJEU ENDIF INA1 = IPLIA(I,1) XDEP = XNORM*Q1(INA1,IND) CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul) XVALA(I,IND,1) = XNORM*XFL XVALA(I,IND,4) = XNORM*XDEP FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL * * ------ choc {l{mentaire POINT_PLAN avec amortissement * ELSE IF (ITYP.EQ.2) THEN XRAID = XPALA(I,1) XJEU = XPALA(I,2) XAMO = XPALA(I,3) IPERM = 0 XNORM = 1. IF (XJEU.LT.XZERO) THEN XNORM = -1. XJEU = -1*XJEU ENDIF INA1 = IPLIA(I,1) XDEP = XNORM * Q1(INA1,IND) IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN XVIT = XNORM * Q2(INA1,IND) ELSE IND2 = IND + 1 XDEPM1 = XNORM * Q1(INA1,IND2) XVIT = (XDEP - XDEPM1) / PDTS2 ENDIF XVALA(I,IND,3) = XNORM*XVIT CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul) *+* IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0 IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0 *+* XVALA(I,IND,1) = XNORM*XFL XVALA(I,IND,4) = XNORM*XDEP FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL * * ------ choc {l{mentaire POINT_PLAN_FLUIDE * ELSE IF (ITYP.EQ.3) THEN XINER = XPALA(I,1) XCONV = XPALA(I,2) XVISC = XPALA(I,3) XPCEL = XPALA(I,4) XPCRA = XPALA(I,5) XJEU = XPALA(I,6) INA1 = IPLIA(I,1) XDEP = Q1(INA1,IND) IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN XVIT = Q2(INA1,IND) ELSE IND2 = IND + 1 XDEPM1 = Q1(INA1,IND2) XVIT = (XDEP - XDEPM1) / PDTS2 ENDIF IF (XJEU.GT.0D0) THEN XDH= XJEU - XDEP XNORM = 1. ELSE XDH= XDEP - XJEU XNORM = -1. ENDIF * Calcul de la masse ajout{e XXIN = -XINER / XDH FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN * Calcul de la force de convection CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul) * Calcul de la force de viscosit{ CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul) * Calcul de la force de perte de charge CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE, & iannul) XFL = (XFCO* XNORM ) + XFVI + XFPE XVALA(I,IND,1) = XDEP XVALA(I,IND,2) = XVIT XVALA(I,IND,3) = XXIN XVALA(I,IND,4) = XFCO*XNORM XVALA(I,IND,5) = XFVI XVALA(I,IND,6) = XFPE FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL * * ------ force {l{mentaire de COUPLAGE EN VITESSE * ELSE IF (ITYP.EQ.4) THEN INA1 = IPLIA(I,1) INA2 = IPLIA(I,2) XDEP = Q1(INA2,IND) IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN XVIT = Q2(INA2,IND) ELSE IND2 = IND + 1 XDEPM1 = Q1(INA2,IND2) XVIT = (XDEP - XDEPM1) / PDTS2 ENDIF XCPLGE = XPALA(I,1) CALL DYCPLV(XVIT,XCPLGE,XFL,iannul) XVALA(I,IND,3) = XVIT XVALA(I,IND,1) = XFL XVALA(I,IND,4) = XDEP FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL * * ------ force {l{mentaire de COUPLAGE EN DEPLACEMENT * ELSE IF (ITYP.EQ.5) THEN INA1 = IPLIA(I,1) INA2 = IPLIA(I,2) XDEP = Q1(INA2,IND) XCPLGE = XPALA(I,1) CALL DYCPLD(XDEP,XCPLGE,XFL,iannul) XVALA(I,IND,1) = XFL XVALA(I,IND,4) = XDEP FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL * * ------ force elementaire de type POLYNOMIALE * ELSE IF (ITYP.EQ.6) THEN * WRITE(IOIMP,*)'DEVLFA : IPLIA(',I,'1)=',IPLIA(I,1) * nombre de modes "origine" NMOD = IPALA(I,2) CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,XDT, & NPAS,I,NMOD,FTOTA,IVINIT,iannul) * * ------ choc {l{mentaire ... * * ELSE IF (ITYP.EQ. ) THEN * ....... * ....... * ENDIF * * la suite n'est plus utile car on passe iannul aux s_p de calcul des * forces de liaisons. * si la liaisn était annulée on l'annule * if ( ( icond.eq. 1 ) .and. ( iannul.eq.1)) then * on annulle l'increment de ftotb * do 112 ik = 1,na1 * do 113 jk = 1,4 * ftota (ik,jk) = ftota0(ik,jk) * 113 continue * 112 continue * end if 10 CONTINUE * END DO * END