1 : C DEVFB8 SOURCE CHAT 05/01/12 22:45:47 5004 2 : C DEVFB7 SOURCE PITO1 97/05/15 13:42:00 2237 3 : SUBROUTINE DEVFB8(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB, 4 : & NPLB,IND,PDT,NUML,IERREU) 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : -INC CCOPTIO 8 : -INC CCREEL 9 : *--------------------------------------------------------------------* 10 : * * 11 : * Opérateur DYNE : algorithme de Fu - de Vogelaere * 12 : * ________________________________________________ * 13 : * * 14 : * Calcul sur la base B des forces hydrodynamiques engendrées * 15 : * dans un palier par le film fluide sur l'arbre * 16 : * (liaisons de type PALIER) * 17 : * * 18 : * Paramètres: * 19 : * * 20 : * e ITYP type de la liaison. * 21 : * es FTOTB Forces extérieures totalisées sur la base B. * 22 : * e XPTB Tableau des déplacements des points. * 23 : * e IPALB Renseigne sur la liaison. * 24 : * e IPLIB Tableau contenant les numéros "DYNE" de la liaison. * 25 : * e XPALB Tableau contenant les paramètres de la liaison. * 26 : * es XVALB Tableau contenant les variables internes de liaisons. * 27 : * e NLIAB Nombre de liaisons sur la base B. * 28 : * e NPLB Nombre total de points intervenant dans les liaisons. * 29 : * e IND Indice du pas. * 30 : * e PDT Durée du pas de temps. * 31 : * e NUML Numéro de la liaison. * 32 : * s IERREU Indicateur d'erreur * 33 : * * 34 : * Auteur, date de création: * 35 : * * 36 : * Valérie BOISSON : le 15 mai 1997 : Création * 37 : * * 38 : *--------------------------------------------------------------------* 39 : * 40 : C On fixe localement les nombres maximum de lobes pour un palier 41 : C et de mailles pour un lobe afin de dimensionner les tableaux 42 : C locaux. Ce choix est arbitraire mais doit permettre la resolution 43 : C de la majorite des problemes reels (on trouve rarement des paliers 44 : C de plus de cinq lobes a ce jour). En cas de modification ulterieure, 45 : C il sera obligatoire de mettre a jour les nouvelles valeurs dans la 46 : C procedure RHODELI. 47 : C 48 : INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*) 49 : REAL*8 XPALB(NLIAB,*),XPTB(NPLB,4,*),FTOTB(NPLB,*) 50 : REAL*8 XVALB(NLIAB,4,*) 51 : LOGICAL ARCPAR 52 : * 53 : * ----- Récupération des paramètres généraux 54 : * 55 : VISCDY = XPALB(NUML,1) 56 : XMASVO = XPALB(NUML,2) 57 : PALM = XPALB(NUML,3) 58 : XLONG = XPALB(NUML,4) 59 : VITROT = XPALB(NUML,7) 60 : PRECIS = XPALB(NUML,10) 61 : NLOBES = IPALB(NUML,6) 62 : NPARRE = IPALB(NUML,7) 63 : * 64 : * ----- Initialisations 65 : * 66 : FNLX=0.D0 67 : FNLY=0.D0 68 : * 69 : ******************************************************************************** 70 : * CAS DES PALIERS CYLINDRIQUES OU A LOBES RESOLU AVEC LE MODELE DE RHODE ET LI 71 : C * 72 : ******************************************************************************** 73 : IF (ITYP.EQ.60) THEN 74 : C NUMP designe le numero "local" du point support de la liaison 75 : NUMP = IPLIB(NUML,1) 76 : 77 : C Soit (R,TETA) le decalage initial du centre du coussinet par rapport 78 : C a l'axe de reference de la structure, exprime en coordonnees polaires : 79 : C R et TETA sont stockees respectivement dans XPALB(NUML,8) et 80 : C XPALB(NUML,9) 81 : * 82 : * ----- Calcul du deplacement et de la vitesse dans le repere local du palier 83 : C 84 : C -------------------> ------------------> --------------------> 85 : C DEPLACEMENT_LOCAL(m) = DECALAGE_COUSSINET + DEPLACEMENT_BASE_B(m) 86 : C 87 : X = XPALB(NUML,8)*COS(XPALB(NUML,9))-XPTB(NUMP,IND,3) 88 : Y = XPALB(NUML,8)*SIN(XPALB(NUML,9))+XPTB(NUMP,IND,2) 89 : C --------------> --------------> 90 : C VITESSE_LOCALE = VITESSE_BASE_B 91 : C 92 : C DEPL(t) - DEPL(t-1) 93 : C = --------------------- 94 : C PDT 95 : C 96 : IF (IND.EQ.2 .OR. IND.EQ.3) THEN 97 : VX = (XPTB(NUMP,IND+1,3)-XPTB(NUMP,IND,3))/PDT*2.D0 98 : VY = (XPTB(NUMP,IND,2)-XPTB(NUMP,IND+1,2))/PDT*2.D0 99 : ELSE 100 : VX = (XPTB(NUMP,IND+2,3)-XPTB(NUMP,IND,3))/PDT 101 : VY = (XPTB(NUMP,IND,2)-XPTB(NUMP,IND+2,2))/PDT 102 : ENDIF 103 : * 104 : * ----- Calcul successif des efforts engendrés par chacun des lobes 105 : * 106 : DO 1 I=1,NLOBES 107 : * 108 : * -- Récupération des paramètres propres aux lobes 109 : * 110 : XJEU = XPALB(NUML,11+NPARRE*(I-1)) 111 : RAYLOB = XPALB(NUML,6)+XPALB(NUML,11+NPARRE*(I-1)) 112 : ASYM = -XPALB(NUML,14+NPARRE*(I-1)) 113 : PRECHA = XPALB(NUML,13+NPARRE*(I-1)) 114 : ANGDEB = XPALB(NUML,12+NPARRE*(I-1))-ASYM 115 : AMPLIT = XPALB(NUML,15+NPARRE*(I-1)) 116 : SURREL = XPALB(NUML,16+NPARRE*(I-1)) 117 : NMAIL = IPALB(NUML,7+I) 118 : c DO 10 J = 1, NMAIL 119 : c G(I,J) = XVALB(NUML,IND+1,7+NMAILP+J) 120 : c G(I,J) = 0 121 : c 10 CONTINUE 122 : * 123 : * -- Calcul de la position et la vitesse du centre de la section 124 : * dans le repère local 125 : * 126 : XAD=(X-PRECHA*COS(ASYM))/XJEU 127 : YAD=(Y-PRECHA*SIN(ASYM))/XJEU 128 : VXAD=VX/XJEU 129 : VYAD=VY/XJEU 130 : IERREU=0 131 : ARCPAR = AMPLIT.LT.(2.D0*XPI) 132 : E=XAD*XAD+YAD*YAD 133 : IF(E.GT.1.D0)GOTO 999 134 : * 135 : * -- Calcul des efforts hydrodynamiques engendres par le film fluide 136 : * -- pour le lobe I 137 : * 137 : CALL RHODLI(XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,I,XLONG, 138 : & VISCDY,RAYLOB,XJEU,VITROT,PRECIS,NMAIL,PALM,XMASVO, 139 : & ARCPAR,SURREL) 140 : * 141 : * -- Totalisation des efforts pour le palier dans sa globalite 142 : * 143 : FNLX=FNLX+FX 144 : FNLY=FNLY+FY 145 : * 146 : * -- Stockage du champ de "pression reduite" G 147 : * 148 : ccccccc 149 : c DO 11 J = 1, NMAIL 150 : c XVALB(NUML,IND,7+NMAILP+J) = G(I,J) 151 : c 11 CONTINUE 152 : c NMAILP = NMAILP + NMAIL 153 : ccccccc 154 : 1 CONTINUE 155 : ******************************************************************************** 156 : * * 157 : * Pour d'autres types de paliers fluides non definis a ce jour ... 158 : C * 159 : * * 160 : ******************************************************************************** 161 : C ELSE IF (ITYP.EQ...) THEN 162 : C 163 : ENDIF 164 : * 165 : * ----- Stockage des variables locales a sortir 166 : * 167 : XVALB(NUML,IND,1) = FNLY 168 : XVALB(NUML,IND,2) = -FNLX 169 : XVALB(NUML,IND,3) = Y 170 : XVALB(NUML,IND,4) = -X 171 : XVALB(NUML,IND,5) = VY 172 : XVALB(NUML,IND,6) = -VX 173 : XVALB(NUML,IND,7) = ((FNLY*VY)+(FNLX*VX))*PDT 174 : IF (XPALB(NUML,5).NE.0) THEN 175 : XPALB(NUML,5) = XPALB(NUML,5) + 0.5D0 176 : IF (IND.EQ.2) WRITE (IOIMP,*) NINT(XPALB(NUML,5)-1) 177 : ENDIF 178 : * 179 : * ----- Totalisation des forces sur la base B 180 : * 181 : FTOTB(NUMP,2) = FTOTB(NUMP,2)+FNLY 182 : FTOTB(NUMP,3) = FTOTB(NUMP,3)-FNLX 183 : GOTO 9999 184 : 185 : C --> Cas ou l'arbre touche le coussinet 186 : 999 IERREU=1 187 : WRITE (IOIMP,*)'EXCENTRICITE SUP A 1' 188 : c -- Fin ajout 189 : XVALB(NUML,IND,1) = 0.D0 190 : XVALB(NUML,IND,2) = 0.D0 191 : c -- Fin ajout 192 : * 193 : * ----- Fin !!! 194 : * 195 : 9999 CONTINUE 196 : RETURN 197 : END 198 : * 199 : ********************************************************************** 200 : ********************************************************************** 201 : * * 202 : * - Sous-programme de calcul des efforts engendrés par un lobe, avec * 203 : * le modèle de Rhode et Li ----------- * 204 : * * 205 : * - Remarque : ce bloc pourra etre separe dans le cas ou de nouveaux * 206 : * modeles de calcul de paliers seraient developpes... * 207 : * * 208 : ********************************************************************** 209 : SUBROUTINE RHODLI(XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,NUMLOB, 210 : & XLONG,VISCDY,RAYLOB,XJEU,VITROT,PRECI,NBMAIL,PALM,XMASVO, 211 : & ARCPAR,SURREL) 212 : * 213 : IMPLICIT INTEGER(I-N) 214 : IMPLICIT REAL*8 (A-H,O-Z) 215 : -INC CCOPTIO 216 : REAL*8 XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,XLONG,VISCDY, 217 : & RAYLOB,XJEU,VITROT,PRECI,PALM,XMASVO,SURREL 218 : C 219 : C Les valeurs de MAXLOB et MAXMAI doivent coincider avec celles definies 220 : C dans l'entete de la subroutine DEVFB8. 221 : C 222 : PARAMETER (MAXLOB = 20) 223 : PARAMETER (MAXMAI = 150) 224 : DIMENSION A(MAXMAI),B(MAXMAI),C(MAXMAI),D(MAXMAI),G(MAXLOB,MAXMAI) 225 : LOGICAL ARCPAR 226 : * 227 : * ----- Initialisations 228 : * 229 : CZ=0.1D-09 230 : ADIM=RAYLOB*RAYLOB/(XJEU*XJEU)*VISCDY*RAYLOB*XLONG 231 : PA0=PALM*XJEU*XJEU/(VISCDY*RAYLOB*RAYLOB) 232 : * 233 : IF (ARCPAR) THEN 234 : G(NUMLOB,1)=0.D0 235 : G(NUMLOB,NBMAIL)=0.D0 236 : DTETA=AMPLIT/DBLE(NBMAIL-1) 237 : ELSE 238 : DTETA=AMPLIT/DBLE(NBMAIL) 239 : ENDIF 240 : c Rajout initialisation G 241 : DO 121 ITRUC2 = 1,MAXMAI 242 : DO 122 ITRUC = 1,MAXLOB 243 : G(ITRUC,ITRUC2)= 0.D0 244 : 122 CONTINUE 245 : 121 CONTINUE 246 : c Fin rajout 247 : 248 : RM=XMASVO/VISCDY*RAYLOB*VITROT 249 : EPSI=SQRT(XAD*XAD+YAD*YAD) 250 : DTCRIT=2.D0*(EPSI*(63.333D0*EPSI-38.D0)+41.2D0) 251 : DO 1 I=1,NBMAIL 252 : * 253 : * ----- Calcul de l'epaisseur du film 254 : * 255 : TETA=DBLE(I-1)*DTETA+ANGDEB 256 : COSTE=COS(TETA) 257 : SINTE=SIN(TETA) 258 : H=1.D0+XAD*COSTE+YAD*SINTE 259 : DHDTET=-XAD*SINTE+YAD*COSTE 260 : * 261 : * ----- Calcul du Reynolds et du Taylor 262 : * 263 : RMH=RM*H*XJEU 264 : TL=RMH*SQRT(H*XJEU/RAYLOB) 265 : IF(TL.LT.DTCRIT) RMH=0.D0 266 : * 267 : * ----- Calcul des coefficients de turbulence 268 : * 269 : AHX=0.0136D0*RMH**0.9D0 270 : AHZ=0.0043D0*RMH**0.96D0 271 : AKX=12.D0+AHX 272 : AKZ=12.D0+AHZ 273 : * 274 : * ----- Calcul des coefficients pour les differences finies 275 : * 276 : H2=H*H 277 : H3=H2*H 278 : DX=DHDTET*H2 279 : 280 : A(I)=(24.D0/(45.D0*DTETA*AKX))*(H3/DTETA+(1.5D0-0.45D0*AHX/AKX) 281 : $ *DX) 282 : B(I)=(-16.D0/3.D0*H3)*(1.D0/(5.D0*AKX*DTETA*DTETA) 283 : & +(RAYLOB/XLONG)*(RAYLOB/XLONG)/AKZ) 284 : C(I)=(24.D0/(45.D0*DTETA*AKX))*(H3/DTETA-(1.5D0-0.45D0*AHX/AKX) 285 : $ *DX) 286 : D(I)=-1.D0/3.D0*(DHDTET*VITROT+2.D0*(VXAD*COSTE+VYAD*SINTE)) 287 : 1 CONTINUE 288 : * 289 : * ----- Calcul par Gauss Seidel 290 : * 291 : 292 : DO 3 ITER=1,500 293 : ER=0.D0 294 : GMAX=0.D0 295 : IF (ARCPAR) THEN 296 : * 297 : * -- Cas du palier partiel (1 lobe avec pression aux deux bords 298 : * materialisant les rainures) 299 : * 300 : DO 4 I=2,NBMAIL-1 301 : S=G(NUMLOB,I) 302 : G(NUMLOB,I)=(1.D0-SURREL)*G(NUMLOB,I)-(SURREL* 303 : & (A(I)*G(NUMLOB,I+1)+C(I)*G(NUMLOB,I-1)+D(I))/B(I)) 304 : IF((G(NUMLOB,I)+PA0).LT.CZ) G(NUMLOB,I)=-PA0 305 : GMAX=MAX(GMAX,G(NUMLOB,I)) 306 : ER=ER+ABS(S-G(NUMLOB,I)) 307 : 4 CONTINUE 308 : ELSE 309 : * 310 : * -- Cas du palier complet (pas de rainure) 311 : * 312 : S=G(NUMLOB,1) 313 : G(NUMLOB,1)=(1.D0-SURREL)*G(NUMLOB,1)-(SURREL* 314 : & (A(1)*G(NUMLOB,2)+C(1)*G(NUMLOB,NBMAIL)+D(1))/B(1)) 315 : IF((G(NUMLOB,1)+PA0).LT.CZ) G(NUMLOB,1)=-PA0 316 : GMAX=MAX(GMAX,G(NUMLOB,1)) 317 : ER=ER+ABS(S-G(NUMLOB,1)) 318 : 319 : S=G(NUMLOB,NBMAIL) 320 : G(NUMLOB,NBMAIL)=(1.D0-SURREL)*G(NUMLOB,NBMAIL)-(SURREL* 321 : & (A(NBMAIL)*G(NUMLOB,1)+C(NBMAIL)*G(NUMLOB,NBMAIL-1) 322 : & +D(NBMAIL))/B(NBMAIL)) 323 : IF((G(NUMLOB,NBMAIL)+PA0).LT.CZ) G(NUMLOB,NBMAIL)=-PA0 324 : GMAX=MAX(GMAX,G(NUMLOB,NBMAIL)) 325 : ER=ER+ABS(S-G(NUMLOB,NBMAIL)) 326 : 327 : DO 44 I=2,NBMAIL-1 328 : S=G(NUMLOB,I) 329 : G(NUMLOB,I)=(1.D0-SURREL)*G(NUMLOB,I)-(SURREL*(A(I)* 330 : & G(NUMLOB,I+1)+C(I)*G(NUMLOB,I-1)+D(I))/B(I)) 331 : IF((G(NUMLOB,I)+PA0).LT.CZ) G(NUMLOB,I)=-PA0 332 : GMAX=MAX(GMAX,G(NUMLOB,I)) 333 : ER=ER+ABS(S-G(NUMLOB,I)) 334 : 335 : 44 CONTINUE 336 : 337 : ENDIF 338 : ERMOY=0.D0 339 : IF (GMAX.GT.CZ) ERMOY=ER/GMAX/DBLE(NBMAIL-2) 340 : IF(ERMOY.LT.PRECI) GOTO 10 341 : 3 CONTINUE 342 : WRITE (IOIMP,*)'attention probleme de convergence! ',ERMOY,GMAX,ER 343 : * 344 : * ----- Calcul des efforts engendres par le lobe (valeur "reduite") 345 : * 346 : 10 FXAD=0.D0 347 : FYAD=0.D0 348 : DO 15 I=1,NBMAIL 349 : TETA=DBLE(I-1)*DTETA+ANGDEB 350 : FXAD=FXAD+G(NUMLOB,I)*COS(TETA) 351 : FYAD=FYAD+G(NUMLOB,I)*SIN(TETA) 352 : 15 CONTINUE 353 : * 354 : * ----- Calcul des efforts engendres par le lobe (valeur dimensionnee) 355 : * 356 : c calcul initialement effectue 357 : FX=FXAD*DTETA*2.D0/3.D0*ADIM 358 : FY=FYAD*DTETA*2.D0/3.D0*ADIM 359 : c a prendre en compte si nous avons un seul patin completement isole 360 : c FX=FXAD*DTETA*2.D0/3.D0*ADIM+PALM*XLONG*RAYLOB*(SIN(AMPLIT+ANGDEB) 361 : c * -SIN(ANGDEB)) 362 : c FY=FYAD*DTETA*2.D0/3.D0*ADIM+PALM*XLONG*RAYLOB*(COS(ANGDEB) 363 : c * -COS(AMPLIT+ANGDEB)) 364 : * 365 : * ----- Fin !!! 366 : * 367 : RETURN 368 : END 369 : 370 : 371 : 372 : 373 :
© Cast3M 2003 - All rights reserved.
Disclaimer