Sources Esope | Cast3M

rhodli
C RHODLI    SOURCE    BP208322  18/07/11    21:15:25     9879           
**********************************************************************
*                                                                    *
* - Sous-programme de calcul des efforts engendrés par un lobe, avec *
* le modèle de Rhode et Li                         -----------       *
* Valérie BOISSON : le 15 mai 1997                                   *
*                                                                    *
**********************************************************************
      SUBROUTINE RHODLI(XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,NUMLOB,
     &     XLONG,VISCDY,RAYLOB,XJEU,VITROT,PRECI,NBMAIL,PALM,XMASVO,
     &     ARCPAR,SURREL,XPROG)
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
      REAL*8 XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,XLONG,VISCDY,
     &       RAYLOB,XJEU,VITROT,PRECI,PALM,XMASVO,SURREL
      REAL*8 XPROG(2,NBMAIL)
C
cbp       PARAMETER (MAXLOB = 20)
cbp       PARAMETER (MAXMAI = 150)
cbp2  PARAMETER (MAXLOB = 10)
cbp2  DIMENSION A(MAXMAI),B(MAXMAI),C(MAXMAI),D(MAXMAI),G(MAXLOB,MAXMAI)
      PARAMETER (MAXMAI = 360)
      DIMENSION A(MAXMAI),B(MAXMAI),C(MAXMAI),D(MAXMAI),G(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
cbp2     G(1)=0.D0
cbp2     G(NBMAIL)=0.D0
         DTETA=AMPLIT/DBLE(NBMAIL-1)
      ELSE
         DTETA=AMPLIT/DBLE(NBMAIL)
      ENDIF
 
c Rajout initialisation G
      DO 121 ITRUC2 = 1,MAXMAI
cbp2         DO 122 ITRUC = 1,MAXLOB
cbp2            G(ITRUC,ITRUC2)= 0.D0
            G(ITRUC2)= 0.D0
cbp2 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
*
cbp2018          TETA=DBLE(I-1)*DTETA+ANGDEB
cbp2018          COSTE=COS(TETA)
cbp2018          SINTE=SIN(TETA)
         COSTE=XPROG(1,I)
         SINTE=XPROG(2,I)
         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(I)
               G(I)=(1.D0-SURREL)*G(I)-(SURREL*
     &              (A(I)*G(I+1)+C(I)*G(I-1)+D(I))/B(I))
               IF((G(I)+PA0).LT.CZ) G(I)=-PA0
               GMAX=MAX(GMAX,G(I))
               ER=ER+ABS(S-G(I))
 4          CONTINUE
 
         ELSE
*
* -- Cas du palier complet (pas de rainure)
*
            S=G(1)
            G(1)=(1.D0-SURREL)*G(1)-(SURREL*
     &           (A(1)*G(2)+C(1)*G(NBMAIL)+D(1))/B(1))
            IF((G(1)+PA0).LT.CZ) G(1)=-PA0
            GMAX=MAX(GMAX,G(1))
            ER=ER+ABS(S-G(1))
 
            S=G(NBMAIL)
            G(NBMAIL)=(1.D0-SURREL)*G(NBMAIL)-(SURREL*
     &           (A(NBMAIL)*G(1)+C(NBMAIL)*G(NBMAIL-1)
     &           +D(NBMAIL))/B(NBMAIL))
            IF((G(NBMAIL)+PA0).LT.CZ) G(NBMAIL)=-PA0
            GMAX=MAX(GMAX,G(NBMAIL))
            ER=ER+ABS(S-G(NBMAIL))
 
            DO 44 I=2,NBMAIL-1
               S=G(I)
               G(I)=(1.D0-SURREL)*G(I)-(SURREL*(A(I)*
     &              G(I+1)+C(I)*G(I-1)+D(I))/B(I))
               IF((G(I)+PA0).LT.CZ) G(I)=-PA0
               GMAX=MAX(GMAX,G(I))
               ER=ER+ABS(S-G(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   CONTINUE
      FXAD=0.D0
      FYAD=0.D0
      DO 15 I=1,NBMAIL
         FXAD=FXAD+(G(I)*XPROG(1,I))
         FYAD=FYAD+(G(I)*XPROG(2,I))
 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