Télécharger rayt1.eso

Retour à la liste

Numérotation des lignes :

  1. C RAYT1 SOURCE CHAT 11/03/16 21:30:26 6902
  2. C RAYE3 SOURCE CHAT 05/01/13 02:45:13 5004
  3. SUBROUTINE RAYT1(MATR, EMIS, TEMP, ERRJ , TRAD, KABS, TABS)
  4.  
  5. C ************************************************************
  6. C **** SUBROUTINE DE CALCUL DE LA TEMPERATURE TRAD ****
  7. C **** INTERVENANT EN RAYONNEMENT THERMIQUE ****
  8. C **** ****
  9. C **** En entree : MATR matrice des facteurs de forme ****
  10. C **** EMIS valeur de l'emissivite en chaque ****
  11. C **** element de surface ****
  12. C **** TEMP temperature moyenne par element ****
  13. C **** ERRJ erreur relative sur le calcul ****
  14. C **** milieu absorbant (02/2011) ****
  15. C **** KABS si =1 il existe un milieu absorbant**
  16. C **** TABS temperaure du milieu absorbant ****
  17. C **** ****
  18. C **** En sortie : TRAD temperature moyenne par element ****
  19. C **** résultant des transferts radiatifs***
  20. C **** ****
  21. C **** ****
  22. C ************************************************************
  23.  
  24. IMPLICIT INTEGER(I-N)
  25. IMPLICIT REAL*8 (A-H,O-Z)
  26.  
  27. -INC CCREEL
  28. -INC CCOPTIO
  29.  
  30. C **********************************************************
  31. C **** Declaration de la structure des facteurs ****
  32. C **** de forme ****
  33. C **********************************************************
  34.  
  35. SEGMENT IFACFO
  36. INTEGER LFACT(NBEL1)
  37. ENDSEGMENT
  38. SEGMENT LFAC
  39. REAL *8 FACT(NBEL2)
  40. ENDSEGMENT
  41.  
  42. C **********************************************************
  43. C **** Declaration des variables du probleme ****
  44. C **********************************************************
  45.  
  46. POINTEUR MATR.IFACFO
  47. POINTEUR LMATR.LFAC
  48. POINTEUR EMIS.LFAC, TEMP.LFAC, TRAD.LFAC
  49. POINTEUR EMIT.LFAC, JRAD.LFAC, JRAD1.LFAC, DIFF.LFAC
  50. C!!
  51. POINTEUR EGAZ.LFAC, AF.LFAC
  52.  
  53. C constante de Stefan
  54. SIG = 5.67D-8
  55. C nombre d'iterations max pour le calcul de la radiosité
  56. NK = 50
  57.  
  58.  
  59. C **********************************************************
  60. C **** Activation de la matrice des facteurs de forme ****
  61. C **** par l'intermediaire de son pointeur. Sa ****
  62. C **** dimension est Nbel. Le dernier pointeur pointe ****
  63. C **** sur les elements de surface ****
  64. C **********************************************************
  65.  
  66. IF (IIMPI.GE.4) THEN
  67. WRITE(6,*) 'DEBUT DE RAYT1.ESO'
  68. ENDIF
  69.  
  70. SEGACT MATR
  71. NEL = MATR.LFACT(/1) - 1
  72.  
  73. NBEL2 = NEL
  74. SEGACT, EMIS, TEMP
  75. SEGINI EMIT, JRAD, JRAD1, TRAD, DIFF
  76. C!!
  77. SEGINI EGAZ, AF
  78.  
  79. CALL UTINIV(DIFF.FACT,NEL)
  80.  
  81. C write(6,*) ' NEL SIG ', NEL, SIG
  82.  
  83. C **********************************************************
  84. C **** calcul des flux emis par les parois ****
  85. C **********************************************************
  86.  
  87. DO I = 1, NEL
  88. EMIT.FACT(I) = EMIS.FACT(I) * SIG * (TEMP.FACT(I)**4)
  89. JRAD1.FACT(I) = EMIT.FACT(I)
  90. ENDDO
  91.  
  92. IF(IIMPI.EQ.(-1)) THEN
  93. write(6,*) ' flux emis'
  94. CALL UTPRIM(EMIT.FACT,NEL)
  95. ENDIF
  96.  
  97. IF (KABS.EQ.1) THEN
  98. C
  99. C **********************************************************
  100. C **** calcul du flux emis par le milieu absorbant ****
  101. C **********************************************************
  102.  
  103. C ------------------------------------------------------
  104. C calcul des sommes sur j des Fij : tableau AF ****
  105. C ------------------------------------------------------
  106. DO I = 1, NEL
  107.  
  108. LMATR = MATR.LFACT(I)
  109. SEGACT LMATR
  110.  
  111. AF.FACT(I) = 0.D0
  112.  
  113. DO J = 1, NEL
  114. AF.FACT(I) = AF.FACT(I) + LMATR.FACT(J)
  115. ENDDO
  116.  
  117. SEGDES LMATR
  118.  
  119. ENDDO
  120.  
  121. C CALL UTPRIM(AF.FACT,NEL)
  122. C ------------------------------------------------------
  123.  
  124. SIGTG4 = SIG * (TABS**4)
  125. C write(6,*) ' SIGTG4: ', SIGTG4
  126.  
  127. DO I = 1, NEL
  128. EGAZ.FACT(I) = (1.D0 - AF.FACT(I))
  129. EGAZ.FACT(I) = EGAZ.FACT(I)*SIGTG4
  130. ENDDO
  131.  
  132.  
  133. IF(IIMPI.EQ.(-1)) THEN
  134. write(6,*) ' flux emis par le gaz'
  135. CALL UTPRIM(EGAZ.FACT,NEL)
  136. ENDIF
  137.  
  138. ENDIF
  139.  
  140. C **********************************************************
  141. C **** calcul iteratif de la radiosité ****
  142. C **** methode de Gauss-Seidel ****
  143. C **********************************************************
  144.  
  145.  
  146. DO 1 K=1,NK
  147.  
  148. C ------------------------------------------------------
  149. DO I = 1, NEL
  150.  
  151. LMATR = MATR.LFACT(I)
  152. SEGACT LMATR
  153.  
  154.  
  155. IF (KABS.EQ.1) THEN
  156.  
  157. JRAD.FACT(I) = EMIT.FACT(I)
  158. $ + (1.D0-EMIS.FACT(I))*EGAZ.FACT(I)
  159.  
  160. ELSE
  161. JRAD.FACT(I) = EMIT.FACT(I)
  162.  
  163. ENDIF
  164.  
  165.  
  166. RO = 1.D0 - EMIS.FACT(I)
  167.  
  168. C cas des surfaces noires: on saute les boucles
  169.  
  170. IF (RO.GT.1D-3) THEN
  171.  
  172.  
  173. C methode de Jacobi pour mémoire
  174.  
  175. C DO J = 1, NEL
  176. C
  177. C JRAD.FACT(I) = JRAD.FACT(I)
  178. C $ + RO * LMATR.FACT(J) * JRAD1.FACT(J)
  179. C ENDDO
  180.  
  181. C... elements tels que J<I
  182. C -------------------------------------------------
  183. IF (I.GT.1) THEN
  184. DO J = 1, (I-1)
  185.  
  186. JRAD.FACT(I) = JRAD.FACT(I)
  187. $ + RO * LMATR.FACT(J) * JRAD.FACT(J)
  188.  
  189. ENDDO
  190. ENDIF
  191. C -------------------------------------------------
  192.  
  193. C... elements tels que J>=I
  194.  
  195. C -------------------------------------------------
  196. DO J = I, NEL
  197.  
  198. JRAD.FACT(I) = JRAD.FACT(I)
  199. $ + RO * LMATR.FACT(J) * JRAD1.FACT(J)
  200. ENDDO
  201. C -------------------------------------------------
  202. C ENDIF
  203.  
  204. ENDIF
  205.  
  206. C JRAD.FACT(I) = EMIT.FACT(I) + JRAD.FACT(I)
  207.  
  208. C COEF = 1.D0 - ( RO * LMATR.FACT(I))
  209. C JRAD.FACT(I) = JRAD.FACT(I) / COEF
  210.  
  211. SEGDES LMATR
  212.  
  213. ENDDO
  214. C ------------------------------------------------------
  215.  
  216. C Test de convergence sur la radiosite
  217.  
  218. DO I = 1, NEL
  219. IF (JRAD1.FACT(I).GE.1D-2) THEN
  220. DIFF.FACT(I) = 1.D0 - (JRAD.FACT(I)/JRAD1.FACT(I))
  221. ENDIF
  222. ENDDO
  223.  
  224. CALL UTMXV(DIFF.FACT,NEL,VAL)
  225. IF (IIMPI.EQ.(-1)) write(6,*) ' K VAL ',K,VAL
  226.  
  227. IF (VAL.LE.ERRJ) THEN
  228. GOTO 2
  229. ENDIF
  230.  
  231. DO I = 1, NEL
  232. JRAD1.FACT(I) = JRAD.FACT(I)
  233. ENDDO
  234.  
  235. IF(IIMPI.EQ.(-1)) THEN
  236. write(6,*) ' JRAD '
  237. CALL UTPRIM(JRAD.FACT,NEL)
  238. write(6,*) ' JRAD1 '
  239. CALL UTPRIM(JRAD1.FACT,NEL)
  240. ENDIF
  241.  
  242. IF(IIMPI.EQ.(-1)) THEN
  243. write(6,*) ' radiosite '
  244. CALL UTPRIM(JRAD.FACT,NEL)
  245. ENDIF
  246.  
  247. 1 CONTINUE
  248. 2 CONTINUE
  249.  
  250. IF (IIMPI.EQ.1) THEN
  251. write (6,*) ' K VAL ',K,VAL
  252. ENDIF
  253.  
  254.  
  255. C **********************************************************
  256. C **** calcul des eclairements ****
  257. C **********************************************************
  258.  
  259. DO I = 1, NEL
  260.  
  261. LMATR = MATR.LFACT(I)
  262. SEGACT LMATR
  263.  
  264. C! CALL UTPRIM(LMATR.FACT,NEL)
  265.  
  266. JRAD.FACT(I) = 0.D0
  267.  
  268. DO J =1, NEL
  269.  
  270. C! write(6,*) ' i j ', I,J
  271. C! write(6,*) 'Fij Jj ', LMATR.FACT(J), JRAD1.FACT(J)
  272.  
  273. JRAD.FACT(I) = JRAD.FACT(I) + LMATR.FACT(J) * JRAD1.FACT(J)
  274. ENDDO
  275.  
  276.  
  277.  
  278. IF (KABS.EQ.1) THEN
  279.  
  280. C on ajoute le terme du au milieu absorbant
  281. C
  282. JRAD.FACT(I) = JRAD.FACT(I) + EGAZ.FACT(I)
  283.  
  284. C! write(6,*) ' Ei ',JRAD.FACT(I)
  285.  
  286. ENDIF
  287.  
  288. SEGDES LMATR
  289. ENDDO
  290.  
  291. IF(IIMPI.EQ.(-1)) THEN
  292. write(6,*) ' eclairement'
  293. CALL UTPRIM(JRAD.FACT,NEL)
  294. ENDIF
  295.  
  296. C **********************************************************
  297. C **** calcul de la temperature TRAD résultat ****
  298. C **********************************************************
  299. C
  300. C étape inchangée si milieu absorbant
  301.  
  302. DO I = 1, NEL
  303. TRAD.FACT(I) = ((1./SIG) * JRAD.FACT(I))**(0.25)
  304. ENDDO
  305.  
  306. IF(IIMPI.EQ.(-1)) THEN
  307. write(6,*) ' Trad'
  308. CALL UTPRIM(TRAD.FACT,NEL)
  309. ENDIF
  310.  
  311. SEGDES EMIS, TEMP, TRAD
  312. SEGDES MATR
  313.  
  314. SEGSUP EMIT, JRAD, JRAD1, DIFF
  315. C!!
  316. SEGSUP EGAZ, AF
  317.  
  318. RETURN
  319. END
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  

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