C DEVFB8 SOURCE CHAT 05/01/12 22:45:47 5004 C DEVFB7 SOURCE PITO1 97/05/15 13:42:00 2237 SUBROUTINE DEVFB8(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB, & NPLB,IND,PDT,NUML,IERREU) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCOPTIO -INC CCREEL *--------------------------------------------------------------------* * * * Opérateur DYNE : algorithme de Fu - de Vogelaere * * ________________________________________________ * * * * Calcul sur la base B des forces hydrodynamiques engendrées * * dans un palier par le film fluide sur l'arbre * * (liaisons de type PALIER) * * * * 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 PDT Durée du pas de temps. * * e NUML Numéro de la liaison. * * s IERREU Indicateur d'erreur * * * * Auteur, date de création: * * * * Valérie BOISSON : le 15 mai 1997 : Création * * * *--------------------------------------------------------------------* * C On fixe localement les nombres maximum de lobes pour un palier C et de mailles pour un lobe afin de dimensionner les tableaux C locaux. Ce choix est arbitraire mais doit permettre la resolution C de la majorite des problemes reels (on trouve rarement des paliers C de plus de cinq lobes a ce jour). En cas de modification ulterieure, C il sera obligatoire de mettre a jour les nouvelles valeurs dans la C procedure RHODELI. C INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*) REAL*8 XPALB(NLIAB,*),XPTB(NPLB,4,*),FTOTB(NPLB,*) REAL*8 XVALB(NLIAB,4,*) LOGICAL ARCPAR * * ----- Récupération des paramètres généraux * VISCDY = XPALB(NUML,1) XMASVO = XPALB(NUML,2) PALM = XPALB(NUML,3) XLONG = XPALB(NUML,4) VITROT = XPALB(NUML,7) PRECIS = XPALB(NUML,10) NLOBES = IPALB(NUML,6) NPARRE = IPALB(NUML,7) * * ----- Initialisations * FNLX=0.D0 FNLY=0.D0 * ******************************************************************************** * CAS DES PALIERS CYLINDRIQUES OU A LOBES RESOLU AVEC LE MODELE DE RHODE ET LI C * ******************************************************************************** IF (ITYP.EQ.60) THEN C NUMP designe le numero "local" du point support de la liaison NUMP = IPLIB(NUML,1) C Soit (R,TETA) le decalage initial du centre du coussinet par rapport C a l'axe de reference de la structure, exprime en coordonnees polaires : C R et TETA sont stockees respectivement dans XPALB(NUML,8) et C XPALB(NUML,9) * * ----- Calcul du deplacement et de la vitesse dans le repere local du palier C C -------------------> ------------------> --------------------> C DEPLACEMENT_LOCAL(m) = DECALAGE_COUSSINET + DEPLACEMENT_BASE_B(m) C X = XPALB(NUML,8)*COS(XPALB(NUML,9))-XPTB(NUMP,IND,3) Y = XPALB(NUML,8)*SIN(XPALB(NUML,9))+XPTB(NUMP,IND,2) C --------------> --------------> C VITESSE_LOCALE = VITESSE_BASE_B C C DEPL(t) - DEPL(t-1) C = --------------------- C PDT C IF (IND.EQ.2 .OR. IND.EQ.3) THEN VX = (XPTB(NUMP,IND+1,3)-XPTB(NUMP,IND,3))/PDT*2.D0 VY = (XPTB(NUMP,IND,2)-XPTB(NUMP,IND+1,2))/PDT*2.D0 ELSE VX = (XPTB(NUMP,IND+2,3)-XPTB(NUMP,IND,3))/PDT VY = (XPTB(NUMP,IND,2)-XPTB(NUMP,IND+2,2))/PDT ENDIF * * ----- Calcul successif des efforts engendrés par chacun des lobes * DO 1 I=1,NLOBES * * -- Récupération des paramètres propres aux lobes * XJEU = XPALB(NUML,11+NPARRE*(I-1)) RAYLOB = XPALB(NUML,6)+XPALB(NUML,11+NPARRE*(I-1)) ASYM = -XPALB(NUML,14+NPARRE*(I-1)) PRECHA = XPALB(NUML,13+NPARRE*(I-1)) ANGDEB = XPALB(NUML,12+NPARRE*(I-1))-ASYM AMPLIT = XPALB(NUML,15+NPARRE*(I-1)) SURREL = XPALB(NUML,16+NPARRE*(I-1)) NMAIL = IPALB(NUML,7+I) c DO 10 J = 1, NMAIL c G(I,J) = XVALB(NUML,IND+1,7+NMAILP+J) c G(I,J) = 0 c 10 CONTINUE * * -- Calcul de la position et la vitesse du centre de la section * dans le repère local * XAD=(X-PRECHA*COS(ASYM))/XJEU YAD=(Y-PRECHA*SIN(ASYM))/XJEU VXAD=VX/XJEU VYAD=VY/XJEU IERREU=0 ARCPAR = AMPLIT.LT.(2.D0*XPI) E=XAD*XAD+YAD*YAD IF(E.GT.1.D0)GOTO 999 * * -- Calcul des efforts hydrodynamiques engendres par le film fluide * -- pour le lobe I * CALL RHODLI(XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,I,XLONG, & VISCDY,RAYLOB,XJEU,VITROT,PRECIS,NMAIL,PALM,XMASVO, & ARCPAR,SURREL) * * -- Totalisation des efforts pour le palier dans sa globalite * FNLX=FNLX+FX FNLY=FNLY+FY * * -- Stockage du champ de "pression reduite" G * ccccccc c DO 11 J = 1, NMAIL c XVALB(NUML,IND,7+NMAILP+J) = G(I,J) c 11 CONTINUE c NMAILP = NMAILP + NMAIL ccccccc 1 CONTINUE ******************************************************************************** * * * Pour d'autres types de paliers fluides non definis a ce jour ... C * * * ******************************************************************************** C ELSE IF (ITYP.EQ...) THEN C ENDIF * * ----- Stockage des variables locales a sortir * XVALB(NUML,IND,1) = FNLY XVALB(NUML,IND,2) = -FNLX XVALB(NUML,IND,3) = Y XVALB(NUML,IND,4) = -X XVALB(NUML,IND,5) = VY XVALB(NUML,IND,6) = -VX XVALB(NUML,IND,7) = ((FNLY*VY)+(FNLX*VX))*PDT IF (XPALB(NUML,5).NE.0) THEN XPALB(NUML,5) = XPALB(NUML,5) + 0.5D0 IF (IND.EQ.2) WRITE (IOIMP,*) NINT(XPALB(NUML,5)-1) ENDIF * * ----- Totalisation des forces sur la base B * FTOTB(NUMP,2) = FTOTB(NUMP,2)+FNLY FTOTB(NUMP,3) = FTOTB(NUMP,3)-FNLX GOTO 9999 C --> Cas ou l'arbre touche le coussinet 999 IERREU=1 WRITE (IOIMP,*)'EXCENTRICITE SUP A 1' c -- Fin ajout XVALB(NUML,IND,1) = 0.D0 XVALB(NUML,IND,2) = 0.D0 c -- Fin ajout * * ----- Fin !!! * 9999 CONTINUE RETURN END * ********************************************************************** ********************************************************************** * * * - Sous-programme de calcul des efforts engendrés par un lobe, avec * * le modèle de Rhode et Li ----------- * * * * - Remarque : ce bloc pourra etre separe dans le cas ou de nouveaux * * modeles de calcul de paliers seraient developpes... * * * ********************************************************************** SUBROUTINE RHODLI(XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,NUMLOB, & XLONG,VISCDY,RAYLOB,XJEU,VITROT,PRECI,NBMAIL,PALM,XMASVO, & ARCPAR,SURREL) * IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCOPTIO REAL*8 XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,XLONG,VISCDY, & RAYLOB,XJEU,VITROT,PRECI,PALM,XMASVO,SURREL C C Les valeurs de MAXLOB et MAXMAI doivent coincider avec celles definies C dans l'entete de la subroutine DEVFB8. C PARAMETER (MAXLOB = 20) PARAMETER (MAXMAI = 150) DIMENSION A(MAXMAI),B(MAXMAI),C(MAXMAI),D(MAXMAI),G(MAXLOB,MAXMAI) LOGICAL ARCPAR * * ----- Initialisations * CZ=0.1D-09 ADIM=RAYLOB*RAYLOB/(XJEU*XJEU)*VISCDY*RAYLOB*XLONG PA0=PALM*XJEU*XJEU/(VISCDY*RAYLOB*RAYLOB) * IF (ARCPAR) THEN G(NUMLOB,1)=0.D0 G(NUMLOB,NBMAIL)=0.D0 DTETA=AMPLIT/DBLE(NBMAIL-1) ELSE DTETA=AMPLIT/DBLE(NBMAIL) ENDIF c Rajout initialisation G DO 121 ITRUC2 = 1,MAXMAI DO 122 ITRUC = 1,MAXLOB G(ITRUC,ITRUC2)= 0.D0 122 CONTINUE 121 CONTINUE c Fin rajout RM=XMASVO/VISCDY*RAYLOB*VITROT EPSI=SQRT(XAD*XAD+YAD*YAD) DTCRIT=2.D0*(EPSI*(63.333D0*EPSI-38.D0)+41.2D0) DO 1 I=1,NBMAIL * * ----- Calcul de l'epaisseur du film * TETA=DBLE(I-1)*DTETA+ANGDEB COSTE=COS(TETA) SINTE=SIN(TETA) H=1.D0+XAD*COSTE+YAD*SINTE DHDTET=-XAD*SINTE+YAD*COSTE * * ----- Calcul du Reynolds et du Taylor * RMH=RM*H*XJEU TL=RMH*SQRT(H*XJEU/RAYLOB) IF(TL.LT.DTCRIT) RMH=0.D0 * * ----- Calcul des coefficients de turbulence * AHX=0.0136D0*RMH**0.9D0 AHZ=0.0043D0*RMH**0.96D0 AKX=12.D0+AHX AKZ=12.D0+AHZ * * ----- Calcul des coefficients pour les differences finies * H2=H*H H3=H2*H DX=DHDTET*H2 A(I)=(24.D0/(45.D0*DTETA*AKX))*(H3/DTETA+(1.5D0-0.45D0*AHX/AKX) $ *DX) B(I)=(-16.D0/3.D0*H3)*(1.D0/(5.D0*AKX*DTETA*DTETA) & +(RAYLOB/XLONG)*(RAYLOB/XLONG)/AKZ) C(I)=(24.D0/(45.D0*DTETA*AKX))*(H3/DTETA-(1.5D0-0.45D0*AHX/AKX) $ *DX) D(I)=-1.D0/3.D0*(DHDTET*VITROT+2.D0*(VXAD*COSTE+VYAD*SINTE)) 1 CONTINUE * * ----- Calcul par Gauss Seidel * DO 3 ITER=1,500 ER=0.D0 GMAX=0.D0 IF (ARCPAR) THEN * * -- Cas du palier partiel (1 lobe avec pression aux deux bords * materialisant les rainures) * DO 4 I=2,NBMAIL-1 S=G(NUMLOB,I) G(NUMLOB,I)=(1.D0-SURREL)*G(NUMLOB,I)-(SURREL* & (A(I)*G(NUMLOB,I+1)+C(I)*G(NUMLOB,I-1)+D(I))/B(I)) IF((G(NUMLOB,I)+PA0).LT.CZ) G(NUMLOB,I)=-PA0 GMAX=MAX(GMAX,G(NUMLOB,I)) ER=ER+ABS(S-G(NUMLOB,I)) 4 CONTINUE ELSE * * -- Cas du palier complet (pas de rainure) * S=G(NUMLOB,1) G(NUMLOB,1)=(1.D0-SURREL)*G(NUMLOB,1)-(SURREL* & (A(1)*G(NUMLOB,2)+C(1)*G(NUMLOB,NBMAIL)+D(1))/B(1)) IF((G(NUMLOB,1)+PA0).LT.CZ) G(NUMLOB,1)=-PA0 GMAX=MAX(GMAX,G(NUMLOB,1)) ER=ER+ABS(S-G(NUMLOB,1)) S=G(NUMLOB,NBMAIL) G(NUMLOB,NBMAIL)=(1.D0-SURREL)*G(NUMLOB,NBMAIL)-(SURREL* & (A(NBMAIL)*G(NUMLOB,1)+C(NBMAIL)*G(NUMLOB,NBMAIL-1) & +D(NBMAIL))/B(NBMAIL)) IF((G(NUMLOB,NBMAIL)+PA0).LT.CZ) G(NUMLOB,NBMAIL)=-PA0 GMAX=MAX(GMAX,G(NUMLOB,NBMAIL)) ER=ER+ABS(S-G(NUMLOB,NBMAIL)) DO 44 I=2,NBMAIL-1 S=G(NUMLOB,I) G(NUMLOB,I)=(1.D0-SURREL)*G(NUMLOB,I)-(SURREL*(A(I)* & G(NUMLOB,I+1)+C(I)*G(NUMLOB,I-1)+D(I))/B(I)) IF((G(NUMLOB,I)+PA0).LT.CZ) G(NUMLOB,I)=-PA0 GMAX=MAX(GMAX,G(NUMLOB,I)) ER=ER+ABS(S-G(NUMLOB,I)) 44 CONTINUE ENDIF ERMOY=0.D0 IF (GMAX.GT.CZ) ERMOY=ER/GMAX/DBLE(NBMAIL-2) IF(ERMOY.LT.PRECI) GOTO 10 3 CONTINUE WRITE (IOIMP,*)'attention probleme de convergence! ',ERMOY,GMAX,ER * * ----- Calcul des efforts engendres par le lobe (valeur "reduite") * 10 FXAD=0.D0 FYAD=0.D0 DO 15 I=1,NBMAIL TETA=DBLE(I-1)*DTETA+ANGDEB FXAD=FXAD+G(NUMLOB,I)*COS(TETA) FYAD=FYAD+G(NUMLOB,I)*SIN(TETA) 15 CONTINUE * * ----- Calcul des efforts engendres par le lobe (valeur dimensionnee) * c calcul initialement effectue FX=FXAD*DTETA*2.D0/3.D0*ADIM FY=FYAD*DTETA*2.D0/3.D0*ADIM c a prendre en compte si nous avons un seul patin completement isole c FX=FXAD*DTETA*2.D0/3.D0*ADIM+PALM*XLONG*RAYLOB*(SIN(AMPLIT+ANGDEB) c * -SIN(ANGDEB)) c FY=FYAD*DTETA*2.D0/3.D0*ADIM+PALM*XLONG*RAYLOB*(COS(ANGDEB) c * -COS(AMPLIT+ANGDEB)) * * ----- Fin !!! * RETURN END