Télécharger rhodli.eso

Retour à la liste

Numérotation des lignes :

  1. C RHODLI SOURCE BP208322 15/07/22 21:15:38 8586
  2. **********************************************************************
  3. * *
  4. * - Sous-programme de calcul des efforts engendrés par un lobe, avec *
  5. * le modèle de Rhode et Li ----------- *
  6. * Valérie BOISSON : le 15 mai 1997 *
  7. * *
  8. **********************************************************************
  9. SUBROUTINE RHODLI(XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,NUMLOB,
  10. & XLONG,VISCDY,RAYLOB,XJEU,VITROT,PRECI,NBMAIL,PALM,XMASVO,
  11. & ARCPAR,SURREL)
  12. *
  13. IMPLICIT INTEGER(I-N)
  14. IMPLICIT REAL*8 (A-H,O-Z)
  15. -INC CCOPTIO
  16. REAL*8 XAD,YAD,VXAD,VYAD,FX,FY,AMPLIT,ANGDEB,XLONG,VISCDY,
  17. & RAYLOB,XJEU,VITROT,PRECI,PALM,XMASVO,SURREL
  18. C
  19. cbp PARAMETER (MAXLOB = 20)
  20. cbp PARAMETER (MAXMAI = 150)
  21. PARAMETER (MAXLOB = 10)
  22. PARAMETER (MAXMAI = 360)
  23. DIMENSION A(MAXMAI),B(MAXMAI),C(MAXMAI),D(MAXMAI),G(MAXLOB,MAXMAI)
  24. LOGICAL ARCPAR
  25. *
  26. * ----- Initialisations
  27. *
  28. CZ=0.1D-09
  29. ADIM=RAYLOB*RAYLOB/(XJEU*XJEU)*VISCDY*RAYLOB*XLONG
  30. PA0=PALM*XJEU*XJEU/(VISCDY*RAYLOB*RAYLOB)
  31. *
  32. IF (ARCPAR) THEN
  33. G(NUMLOB,1)=0.D0
  34. G(NUMLOB,NBMAIL)=0.D0
  35. DTETA=AMPLIT/DBLE(NBMAIL-1)
  36. ELSE
  37. DTETA=AMPLIT/DBLE(NBMAIL)
  38. ENDIF
  39. c Rajout initialisation G
  40. DO 121 ITRUC2 = 1,MAXMAI
  41. DO 122 ITRUC = 1,MAXLOB
  42. G(ITRUC,ITRUC2)= 0.D0
  43. 122 CONTINUE
  44. 121 CONTINUE
  45. c Fin rajout
  46.  
  47. RM=XMASVO/VISCDY*RAYLOB*VITROT
  48. EPSI=SQRT(XAD*XAD+YAD*YAD)
  49. DTCRIT=2.D0*(EPSI*(63.333D0*EPSI-38.D0)+41.2D0)
  50. DO 1 I=1,NBMAIL
  51. *
  52. * ----- Calcul de l'epaisseur du film
  53. *
  54. TETA=DBLE(I-1)*DTETA+ANGDEB
  55. COSTE=COS(TETA)
  56. SINTE=SIN(TETA)
  57. H=1.D0+XAD*COSTE+YAD*SINTE
  58. DHDTET=-XAD*SINTE+YAD*COSTE
  59. *
  60. * ----- Calcul du Reynolds et du Taylor
  61. *
  62. RMH=RM*H*XJEU
  63. TL=RMH*SQRT(H*XJEU/RAYLOB)
  64. IF(TL.LT.DTCRIT) RMH=0.D0
  65. *
  66. * ----- Calcul des coefficients de turbulence
  67. *
  68. AHX=0.0136D0*RMH**0.9D0
  69. AHZ=0.0043D0*RMH**0.96D0
  70. AKX=12.D0+AHX
  71. AKZ=12.D0+AHZ
  72. *
  73. * ----- Calcul des coefficients pour les differences finies
  74. *
  75. H2=H*H
  76. H3=H2*H
  77. DX=DHDTET*H2
  78.  
  79. A(I)=(24.D0/(45.D0*DTETA*AKX))*(H3/DTETA+(1.5D0-0.45D0*AHX/AKX)
  80. $ *DX)
  81. B(I)=(-16.D0/3.D0*H3)*(1.D0/(5.D0*AKX*DTETA*DTETA)
  82. & +(RAYLOB/XLONG)*(RAYLOB/XLONG)/AKZ)
  83. C(I)=(24.D0/(45.D0*DTETA*AKX))*(H3/DTETA-(1.5D0-0.45D0*AHX/AKX)
  84. $ *DX)
  85. D(I)=-1.D0/3.D0*(DHDTET*VITROT+2.D0*(VXAD*COSTE+VYAD*SINTE))
  86. 1 CONTINUE
  87. *
  88. * ----- Calcul par Gauss Seidel
  89. *
  90.  
  91. DO 3 ITER=1,500
  92. ER=0.D0
  93. GMAX=0.D0
  94. IF (ARCPAR) THEN
  95. *
  96. * -- Cas du palier partiel (1 lobe avec pression aux deux bords
  97. * materialisant les rainures)
  98. *
  99. DO 4 I=2,NBMAIL-1
  100. S=G(NUMLOB,I)
  101. G(NUMLOB,I)=(1.D0-SURREL)*G(NUMLOB,I)-(SURREL*
  102. & (A(I)*G(NUMLOB,I+1)+C(I)*G(NUMLOB,I-1)+D(I))/B(I))
  103. IF((G(NUMLOB,I)+PA0).LT.CZ) G(NUMLOB,I)=-PA0
  104. GMAX=MAX(GMAX,G(NUMLOB,I))
  105. ER=ER+ABS(S-G(NUMLOB,I))
  106. 4 CONTINUE
  107. ELSE
  108. *
  109. * -- Cas du palier complet (pas de rainure)
  110. *
  111. S=G(NUMLOB,1)
  112. G(NUMLOB,1)=(1.D0-SURREL)*G(NUMLOB,1)-(SURREL*
  113. & (A(1)*G(NUMLOB,2)+C(1)*G(NUMLOB,NBMAIL)+D(1))/B(1))
  114. IF((G(NUMLOB,1)+PA0).LT.CZ) G(NUMLOB,1)=-PA0
  115. GMAX=MAX(GMAX,G(NUMLOB,1))
  116. ER=ER+ABS(S-G(NUMLOB,1))
  117.  
  118. S=G(NUMLOB,NBMAIL)
  119. G(NUMLOB,NBMAIL)=(1.D0-SURREL)*G(NUMLOB,NBMAIL)-(SURREL*
  120. & (A(NBMAIL)*G(NUMLOB,1)+C(NBMAIL)*G(NUMLOB,NBMAIL-1)
  121. & +D(NBMAIL))/B(NBMAIL))
  122. IF((G(NUMLOB,NBMAIL)+PA0).LT.CZ) G(NUMLOB,NBMAIL)=-PA0
  123. GMAX=MAX(GMAX,G(NUMLOB,NBMAIL))
  124. ER=ER+ABS(S-G(NUMLOB,NBMAIL))
  125.  
  126. DO 44 I=2,NBMAIL-1
  127. S=G(NUMLOB,I)
  128. G(NUMLOB,I)=(1.D0-SURREL)*G(NUMLOB,I)-(SURREL*(A(I)*
  129. & G(NUMLOB,I+1)+C(I)*G(NUMLOB,I-1)+D(I))/B(I))
  130. IF((G(NUMLOB,I)+PA0).LT.CZ) G(NUMLOB,I)=-PA0
  131. GMAX=MAX(GMAX,G(NUMLOB,I))
  132. ER=ER+ABS(S-G(NUMLOB,I))
  133.  
  134. 44 CONTINUE
  135.  
  136. ENDIF
  137. ERMOY=0.D0
  138. IF (GMAX.GT.CZ) ERMOY=ER/GMAX/DBLE(NBMAIL-2)
  139. IF(ERMOY.LT.PRECI) GOTO 10
  140. 3 CONTINUE
  141. WRITE (IOIMP,*)'attention probleme de convergence! ',ERMOY,GMAX,ER
  142. *
  143. * ----- Calcul des efforts engendres par le lobe (valeur "reduite")
  144. *
  145. 10 FXAD=0.D0
  146. FYAD=0.D0
  147. DO 15 I=1,NBMAIL
  148. TETA=DBLE(I-1)*DTETA+ANGDEB
  149. FXAD=FXAD+G(NUMLOB,I)*COS(TETA)
  150. FYAD=FYAD+G(NUMLOB,I)*SIN(TETA)
  151. 15 CONTINUE
  152. *
  153. * ----- Calcul des efforts engendres par le lobe (valeur dimensionnee)
  154. *
  155. c calcul initialement effectue
  156. FX=FXAD*DTETA*2.D0/3.D0*ADIM
  157. FY=FYAD*DTETA*2.D0/3.D0*ADIM
  158. c a prendre en compte si nous avons un seul patin completement isole
  159. c FX=FXAD*DTETA*2.D0/3.D0*ADIM+PALM*XLONG*RAYLOB*(SIN(AMPLIT+ANGDEB)
  160. c * -SIN(ANGDEB))
  161. c FY=FYAD*DTETA*2.D0/3.D0*ADIM+PALM*XLONG*RAYLOB*(COS(ANGDEB)
  162. c * -COS(AMPLIT+ANGDEB))
  163. *
  164. * ----- Fin !!!
  165. *
  166. RETURN
  167. END
  168.  
  169.  
  170.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales