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 : 78 : CALL UTINIV(DIFF.FACT,NEL) 79 : 80 : C write(6,*) ' NEL SIG ', NEL, SIG 81 : 82 : C ********************************************************** 83 : C **** calcul des flux emis par les parois **** 84 : C ********************************************************** 85 : 86 : DO I = 1, NEL 87 : EMIT.FACT(I) = EMIS.FACT(I) * SIG * (TEMP.FACT(I)**4) 88 : JRAD1.FACT(I) = EMIT.FACT(I) 89 : ENDDO 90 : 91 : IF(IIMPI.EQ.(-1)) THEN 92 : write(6,*) ' flux emis' 92 : CALL UTPRIM(EMIT.FACT,NEL) 93 : ENDIF 94 : 95 : IF (KABS.EQ.1) THEN 96 : C 97 : C ********************************************************** 98 : C **** calcul du flux emis par le milieu absorbant **** 99 : C ********************************************************** 100 : 101 : C ------------------------------------------------------ 102 : C calcul des sommes sur j des Fij : tableau AF **** 103 : C ------------------------------------------------------ 104 : DO I = 1, NEL 105 : 106 : LMATR = MATR.LFACT(I) 107 : SEGACT LMATR 108 : 109 : AF.FACT(I) = 0.D0 110 : 111 : DO J = 1, NEL 112 : AF.FACT(I) = AF.FACT(I) + LMATR.FACT(J) 113 : ENDDO 114 : 115 : SEGDES LMATR 116 : 117 : ENDDO 118 : 118 : C CALL UTPRIM(AF.FACT,NEL) 119 : C ------------------------------------------------------ 120 : 121 : SIGTG4 = SIG * (TABS**4) 122 : C write(6,*) ' SIGTG4: ', SIGTG4 123 : 124 : DO I = 1, NEL 125 : EGAZ.FACT(I) = (1.D0 - AF.FACT(I)) 126 : EGAZ.FACT(I) = EGAZ.FACT(I)*SIGTG4 127 : ENDDO 128 : 129 : 130 : IF(IIMPI.EQ.(-1)) THEN 131 : write(6,*) ' flux emis par le gaz' 131 : CALL UTPRIM(EGAZ.FACT,NEL) 132 : ENDIF 133 : 134 : ENDIF 135 : 136 : C ********************************************************** 137 : C **** calcul iteratif de la radiosité **** 138 : C **** methode de Gauss-Seidel **** 139 : C ********************************************************** 140 : 141 : 142 : DO 1 K=1,NK 143 : 144 : C ------------------------------------------------------ 145 : DO I = 1, NEL 146 : 147 : LMATR = MATR.LFACT(I) 148 : SEGACT LMATR 149 : 150 : 151 : IF (KABS.EQ.1) THEN 152 : 153 : JRAD.FACT(I) = EMIT.FACT(I) 154 : $ + (1.D0-EMIS.FACT(I))*EGAZ.FACT(I) 155 : 156 : ELSE 157 : JRAD.FACT(I) = EMIT.FACT(I) 158 : 159 : ENDIF 160 : 161 : 162 : RO = 1.D0 - EMIS.FACT(I) 163 : 164 : C cas des surfaces noires: on saute les boucles 165 : 166 : IF (RO.GT.1D-3) THEN 167 : 168 : 169 : C methode de Jacobi pour mémoire 170 : 171 : C DO J = 1, NEL 172 : C 173 : C JRAD.FACT(I) = JRAD.FACT(I) 174 : C $ + RO * LMATR.FACT(J) * JRAD1.FACT(J) 175 : C ENDDO 176 : 177 : C... elements tels que J=I 190 : 191 : C ------------------------------------------------- 192 : DO J = I, NEL 193 : 194 : JRAD.FACT(I) = JRAD.FACT(I) 195 : $ + RO * LMATR.FACT(J) * JRAD1.FACT(J) 196 : ENDDO 197 : C ------------------------------------------------- 198 : C ENDIF 199 : 200 : ENDIF 201 : 202 : C JRAD.FACT(I) = EMIT.FACT(I) + JRAD.FACT(I) 203 : 204 : C COEF = 1.D0 - ( RO * LMATR.FACT(I)) 205 : C JRAD.FACT(I) = JRAD.FACT(I) / COEF 206 : 207 : SEGDES LMATR 208 : 209 : ENDDO 210 : C ------------------------------------------------------ 211 : 212 : C Test de convergence sur la radiosite 213 : 214 : DO I = 1, NEL 215 : IF (JRAD1.FACT(I).GE.1D-2) THEN 216 : DIFF.FACT(I) = 1.D0 - (JRAD.FACT(I)/JRAD1.FACT(I)) 217 : ENDIF 218 : ENDDO 219 : 219 : CALL UTMXV(DIFF.FACT,NEL,VAL) 220 : IF (IIMPI.EQ.(-1)) write(6,*) ' K VAL ',K,VAL 221 : 222 : IF (VAL.LE.ERRJ) THEN 223 : GOTO 2 224 : ENDIF 225 : 226 : DO I = 1, NEL 227 : JRAD1.FACT(I) = JRAD.FACT(I) 228 : ENDDO 229 : 230 : IF(IIMPI.EQ.(-1)) THEN 231 : write(6,*) ' JRAD ' 231 : CALL UTPRIM(JRAD.FACT,NEL) 232 : write(6,*) ' JRAD1 ' 232 : CALL UTPRIM(JRAD1.FACT,NEL) 233 : ENDIF 234 : 235 : IF(IIMPI.EQ.(-1)) THEN 236 : write(6,*) ' radiosite ' 236 : CALL UTPRIM(JRAD.FACT,NEL) 237 : ENDIF 238 : 239 : 1 CONTINUE 240 : 2 CONTINUE 241 : 242 : IF (IIMPI.EQ.1) THEN 243 : write (6,*) ' K VAL ',K,VAL 244 : ENDIF 245 : 246 : 247 : C ********************************************************** 248 : C **** calcul des eclairements **** 249 : C ********************************************************** 250 : 251 : DO I = 1, NEL 252 : 253 : LMATR = MATR.LFACT(I) 254 : SEGACT LMATR 255 : 255 : C! CALL UTPRIM(LMATR.FACT,NEL) 256 : 257 : JRAD.FACT(I) = 0.D0 258 : 259 : DO J =1, NEL 260 : 261 : C! write(6,*) ' i j ', I,J 262 : C! write(6,*) 'Fij Jj ', LMATR.FACT(J), JRAD1.FACT(J) 263 : 264 : JRAD.FACT(I) = JRAD.FACT(I) + LMATR.FACT(J) * JRAD1.FACT(J) 265 : ENDDO 266 : 267 : 268 : 269 : IF (KABS.EQ.1) THEN 270 : 271 : C on ajoute le terme du au milieu absorbant 272 : C 273 : JRAD.FACT(I) = JRAD.FACT(I) + EGAZ.FACT(I) 274 : 275 : C! write(6,*) ' Ei ',JRAD.FACT(I) 276 : 277 : ENDIF 278 : 279 : SEGDES LMATR 280 : ENDDO 281 : 282 : IF(IIMPI.EQ.(-1)) THEN 283 : write(6,*) ' eclairement' 283 : CALL UTPRIM(JRAD.FACT,NEL) 284 : ENDIF 285 : 286 : C ********************************************************** 287 : C **** calcul de la temperature TRAD résultat **** 288 : C ********************************************************** 289 : C 290 : C étape inchangée si milieu absorbant 291 : 292 : DO I = 1, NEL 293 : TRAD.FACT(I) = ((1./SIG) * JRAD.FACT(I))**(0.25) 294 : ENDDO 295 : 296 : IF(IIMPI.EQ.(-1)) THEN 297 : write(6,*) ' Trad' 297 : CALL UTPRIM(TRAD.FACT,NEL) 298 : ENDIF 299 : 300 : SEGDES EMIS, TEMP, TRAD 301 : SEGDES MATR 302 : 303 : SEGSUP EMIT, JRAD, JRAD1, DIFF 304 : C!! 305 : SEGSUP EGAZ, AF 306 : 307 : RETURN 308 : END 309 : 310 : 311 : 312 : 313 : 314 : 315 :
© Cast3M 2003 - All rights reserved.
Disclaimer