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

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