C DEVFB4 SOURCE BP208322 20/09/18 21:15:29 10718 SUBROUTINE DEVFB4(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB, & NPLB,IND,IND1,INDM1,NPA,NPAM1,IND2,PDT,PDTS2, & FEXPSM,NPC1,I,IERRD,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 PROFIL_PROFIL_... * * * * 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. * * es IERRD num{ro d'erreur. * * * * * * Auteur, date de cr{ation: * * * * Lionel VIVAN : le 1 f{vrier 1990 * * * *--------------------------------------------------------------------* * INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*) REAL*8 XPALB(NLIAB,*),XPTB(NPLB,2,*),FTOTB(NPLB,*) REAL*8 XVALB(NLIAB,4,*),XFOR(3),FEXPSM(NPLB,NPC1,2,*) LOGICAL LPETPP PARAMETER ( ZERO = 0.D0 , PRECIS = 1.D-15 ) * * --- choc {l{mentaire PROFIL_PROFIL_INTERIEUR * --- choc {l{mentaire PROFIL_PROFIL_EXTERIEUR * IF (ITYP.EQ.31 .OR. ITYP.EQ.32) THEN NPOI = IPLIB(I,1) IDIM = IPALB(I,3) NOMBN1 = IPALB(I,4) NOMBN2 = IPALB(I,5) XRAID = XPALB(I,1) XSECT = XPALB(I,2) XEXPO = XPALB(I,3) ID1 = 3 ID2 = ID1 + IDIM ID3 = ID1 + 2*IDIM ID4 = ID1 + 3*IDIM ID5 = ID1 + 4*IDIM ID6 = ID1 + 5*IDIM ID7 = ID6 + IDIM*NOMBN1 ID8 = ID7 + IDIM*NOMBN2 IP1 = 5 IP2 = IP1 + NOMBN1 IP3 = IP2 IP4 = IP1 + NOMBN1 + NOMBN1*NOMBN2 SOMFN = ZERO * * Recherche du nombre de point qui ont travers{ le maillage fixe, * donc s'il y a eu choc. IPOIN = 0 ITIN = ID7 DO 10 IN = 1,NOMBN2 IP5 = IP4 + NOMBN1*(IN - 1) * D{placement du point IN du profil mobile dans le plan des profils XN1 = ZERO YN1 = ZERO DO 12 ID = 1,IDIM XXX = XPALB(I,ITIN+ID) + XPTB(NPOI,1,ID) & - XPALB(I,ID2+ID) XN1 = XN1 + ( XXX * XPALB(I,ID3+ID) ) YN1 = YN1 + ( XXX * XPALB(I,ID4+ID) ) 12 CONTINUE * end do * On boucle sur le nombre d'{l{ment du profil fixe pour savoir si * le point IN @ travers{ le profil fixe NUME1 = ID6 NUME2 = ID6 + IDIM ICOEF = ID8 DO 14 IE = 1,NOMBN1 IPALB(I,IP5+IE) = 0 IPOSI = IPALB(I,IP3+IE) ICAS = IPALB(I,IP1+IE) XAIJ = XPALB(I,ICOEF+1) XBIJ = XPALB(I,ICOEF+2) IF (ICAS.EQ.1) THEN IF (XN1.GT.XAIJ) THEN IPOS2 = 1 ELSE IF (XN1.LT.XAIJ) THEN IPOS2 = -1 ELSE IPOS2 = 0 ENDIF ELSE IF (ICAS.EQ.2) THEN IF (YN1.GT.XBIJ) THEN IPOS2 = 1 ELSE IF (YN1.LT.XBIJ) THEN IPOS2 = -1 ELSE IPOS2 = 0 ENDIF ELSE RES = ( XAIJ * XN1 ) + XBIJ IF (YN1.GT.RES) THEN IPOS2 = 1 ELSE IF (YN1.LT.RES) THEN IPOS2 = -1 ELSE IPOS2 = 0 ENDIF ENDIF IF (IPOSI.NE.IPOS2) THEN * Le point IN a travers{ le maillage fixe, * quel {l{ment a {t{ travers{ par ce point ? IF (IE.EQ.NOMBN1) NUME2 = ID6 XE1 = ZERO YE1 = ZERO XE2 = ZERO YE2 = ZERO DO 16 ID = 1,IDIM XXX = XPALB(I,NUME1+ID) - XPALB(I,ID2+ID) YYY = XPALB(I,NUME2+ID) - XPALB(I,ID2+ID) XE1 = XE1 + ( XXX * XPALB(I,ID3+ID) ) YE1 = YE1 + ( XXX * XPALB(I,ID4+ID) ) XE2 = XE2 + ( YYY * XPALB(I,ID3+ID) ) YE2 = YE2 + ( YYY * XPALB(I,ID4+ID) ) 16 CONTINUE * end do IF (XE2.GT.XE1) THEN IF (XE1.LE.XN1 .AND. XN1.LT.XE2) THEN IPOIN = IPOIN + 1 IPALB(I,IP5+IE) = 1 GOTO 18 ENDIF ELSE IF (XE2.LT.XN1 .AND. XN1.LE.XE1) THEN IPOIN = IPOIN + 1 IPALB(I,IP5+IE) = 1 GOTO 18 ENDIF ENDIF IF (YE2.GT.YE1) THEN IF (YE1.LE.YN1 .AND. YN1.LT.YE2) THEN IPOIN = IPOIN + 1 IPALB(I,IP5+IE) = 1 GOTO 18 ENDIF ELSE IF (YE2.LT.YN1 .AND. YN1.LE.YE1) THEN IPOIN = IPOIN + 1 IPALB(I,IP5+IE) = 1 GOTO 18 ENDIF ENDIF ENDIF NUME1 = NUME2 NUME2 = NUME2 + IDIM ICOEF = ICOEF + 2 14 CONTINUE * end do 18 CONTINUE IP3 = IP3 + NOMBN1 ITIN = ITIN + IDIM 10 CONTINUE * end do * * Si le nombre de points ayant travers{s le profil fixe est {gale * au nombre de points du maillage mobile, alors on imprime une * erreur et on arrete le calcul. IF (IPOIN.EQ.NOMBN2) THEN IERRD = 1 RETURN ENDIF * * s'il y a eu choc, * il y a IPOIN qui ont travers{ le maillage fixe * IF (IPOIN.NE.0) THEN * * Il y a eu choc : * - recherche du premier point qui a travers{ le maillage fixe, d'o: * d{termination du point P, intersection maillage fixe-mobile * - recherche du dernier point qui a travers{ le maillage fixe, d'o: * d{termination du point PP, intersection maillage fixe-mobile * - @ partir des points P et PP, d{termination de la normale de choc * - calcul de l'aire de choc, @ partir du point P, des points du * maillage fixe, du point PP et des points du maillage mobile qui * n'ont pas travers{ le maillage fixe * - calcul de la force de choc. * ICOMP = 0 IN2 = 0 20 CONTINUE * On boucle sur les points du profil mobile IN2 = IN2 + 1 IP5 = IP4 + NOMBN1*(IN2 - 1) IF (IN2.GT.NOMBN2) GOTO 22 LPETPP = .FALSE. DO 24 IE1 = 1,NOMBN1 IETAT = IPALB(I,IP5+IE1) * On regarde si le point IN2 a travers{ un {l{ment du profil fixe IF (IETAT.EQ.1) THEN ICOMP = ICOMP + 1 IF (IN2.EQ.1) THEN * Si on traite le premier point, on cherche les points en partant du * dernier DO 26 INN2 = NOMBN2,1,-1 IP6 = IP4 + NOMBN1*(INN2 - 1) DO 28 IEE1 = 1,NOMBN1 IETA2 = IPALB(I,IP6+IEE1) IF (IETA2.EQ.1) THEN ICOMP = ICOMP + 1 IEEE1 = IEE1 GOTO 26 ENDIF 28 CONTINUE * end do * Le point INN2 n'a pas travers{ le maillage fixe * Le point (INN2-1) et la point INN2 forment le segment de droite * pour d{finir le point d'intersection P avec l'{l{ment IEEE1 IF (INN2.EQ.NOMBN2) THEN IPOIN1 = NOMBN2 IPOIN2 = 1 IELEM1 = IE1 ELSE IPOIN1 = INN2 IPOIN2 = INN2 + 1 IELEM1 = IEEE1 ENDIF GOTO 30 26 CONTINUE * end do 30 CONTINUE ELSE IPOIN1 = IN2 - 1 IPOIN2 = IN2 IELEM1 = IE1 ENDIF * D{termination du point P * Le point IPOIN1 n'a pas travers{ le maillage fixe * Le point IPOIN2 a travers{ le maillage fixe par l'{l{ment IELEM1 CALL DYNE31(IPOIN1,IPOIN2,IELEM1,XPALB,IPALB,XPTB, & NLIAB,NPLB,I,NPOI,IND,ID1,IP1,XP,YP) * * Recherche du point PP, on regarde les points suivant IF (IN2.EQ.NOMBN2) THEN IPOIN3 = IN2 IPOIN4 = 1 IELEM2 = IE1 GOTO 32 ENDIF IN22 = IN2 + 1 IEEE1 = IE1 DO 34 INN2 = IN22,NOMBN2 IP7 = IP4 + NOMBN1*(INN2 - 1) DO 36 IEE1 = 1,NOMBN1 IETA2 = IPALB(I,IP7+IEE1) IF (IETA2.EQ.1) THEN ICOMP = ICOMP + 1 IEEE1 = IEE1 IF (INN2.EQ.NOMBN2) THEN IPOIN3 = INN2 IPOIN4 = 1 IELEM2 = IEE1 ENDIF GOTO 34 ENDIF 36 CONTINUE * end do * Le point INN2 n'a pas travers{ le maillage fixe * Le point (INN2-1) et la point INN2 forment le segment de droite * pour d{finir le point d'intersection P avec l'{l{ment IEEE1 IPOIN3 = INN2 - 1 IPOIN4 = INN2 IELEM2 = IEEE1 GOTO 32 34 CONTINUE * end do 32 CONTINUE * D{termination du point PP * Le point IPOIN3 a travers{ le maillage fixe par l'{l{ment IELEM2 * Le point IPOIN4 n'a pas travers{ le maillage fixe CALL DYNE31(IPOIN3,IPOIN4,IELEM2,XPALB,IPALB,XPTB, & NLIAB,NPLB,I,NPOI,IND,ID1,IP1,XPP,YPP) LPETPP = .TRUE. GOTO 40 ENDIF 24 CONTINUE * end do 40 CONTINUE * IF ( LPETPP ) THEN * On a d{termin{ les points P et PP, calcul de la normale XNORM1 = ( XP + XPP ) * 0.5D0 YNORM1 = ( YP + YPP ) * 0.5D0 XXX = XP - XPP YYY = YP - YPP IPT3 = ID7 + IDIM*(IPOIN3-1) XN4 = ZERO YN4 = ZERO DO 42 ID = 1,IDIM XXX = XPALB(I,IPT3+ID) + XPTB(NPOI,1,ID) & - XPALB(I,ID2+ID) XN4 = XN4 + ( XXX * XPALB(I,ID3+ID) ) YN4 = YN4 + ( XXX * XPALB(I,ID4+ID) ) 42 CONTINUE * END DO IF (ABS(YYY).LT.PRECIS) THEN YNN = YN4 - YNORM1 XNORM = ZERO YNORM = YNN / ABS(YNN) ELSE IF (ABS(XXX).LT.PRECIS) THEN XNN = XN4 - XNORM1 XNORM = XNN / ABS(XNN) YNORM = ZERO ELSE * le premier vecteur est P-PP * la normale de choc sera le produit vectoriel P_PP, N XPPP = XPP - XP YPPP = YPP - YP PS = SQRT( XPPP*XPPP + YPPP*YPPP ) XNORM = YPPP * XPALB(I,ID5+3) / PS YNORM = -XPPP * XPALB(I,ID5+3) / PS ENDIF * Calcul de l'aire CALL DYNE32(IPOIN1,IPOIN4,IELEM1,IELEM2,XP,YP,XPP,YPP, & XPALB,IPALB,XPTB,NLIAB,NPLB,I,NPOI,ID1,IP1,IND,SOMAIR) AIRE = XSECT - SOMAIR * Calcul de la force de choc * Liaison conditionnelle IF (iannul.EQ.0) THEN FN = -XRAID * ( AIRE ** XEXPO ) ELSE FN = 0 ENDIF SOMFN = SOMFN + FN * On revient dans le repere global XFOR1 = FN * XNORM XFOR2 = FN * YNORM DO 38 ID = 1,IDIM XXX = XFOR1*XPALB(I,ID3+ID) + XFOR2*XPALB(I,ID4+ID) FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XXX 38 CONTINUE * end do IN2 = INN2 ENDIF IF (ICOMP.EQ.IPOIN) GOTO 22 GOTO 20 22 CONTINUE ENDIF DO 100 ID = 1,IDIM XVALB(I,IND,3+ID) = XPTB(NPOI,1,ID) 100 CONTINUE * end do XVALB(I,IND,1) = SOMFN * * --- choc ........... * * ELSE IF (ITYP.EQ. ) THEN * ....... * ....... * ENDIF * END