Télécharger rhodli.eso

Retour à la liste

Numérotation des lignes :

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

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