Télécharger la plaquette Cast3M


fr en

Télécharger rayt1.eso

Retour à la liste

Numérotation des lignes : 

rayt1

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

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