1 : C BSUR2 SOURCE DURA 06/07/11 21:17:22 5481 2 : SUBROUTINE BSUR2 (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL, 3 : & RINDEX,P1,T1,QAE,QEE1,PHI1,P2,T2,U2,QEE2,PHI2,QW2,RE,H, 4 : & PSQ,NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4) 5 : 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8 (A-H,O-Z) 8 : DIMENSION XN(NPP),TN(NPP),EN(NPP),BN(NPP) 9 : 10 : C operateur FUITE 11 : C>>> superheated vapour 12 : C version with wall condensation 13 : C version gaz reel 14 : C QAE: QA.E.B constant air flowrate 15 : C QEEi: QEi.E.B water flowrate 16 : C indice 1,2 : entree et sortie de troncon 17 : C RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur 18 : C 19 : C H correspond au flux d'energie total cede a la paroi 20 : C HM (convection) est utilise pour le calcul de T 21 : C 22 : IF(KIMP.GE.1) THEN 23 : write(6,*) 'entree bsur2' 24 : ENDIF 25 : 26 : HW=0.D0 27 : RINDEX=1.0 28 : DQEE = 0.D0 29 : 30 : C IMPLEMENTATION DES PROPRIETES PHYSIQUES 30 : CALL BPHYS(T0,P0,RA,RV,CA,XCV,XCL,XLAT0,XROL,XKL,XKT,XREL) 31 : 32 : C TEMPERATURE PAROI A ENTREE (KELVIN) 32 : CALL BTPAR(XN,TN,X,TP,NPP,ITP) 33 : TP=TP+T0 34 : 35 : C EPAISSEUR DE FISSURE EN ENTREE 35 : CALL BTPAR(XN,EN,X,E1,NPP,ITP) 36 : C etendue a l'entree 36 : CALL BTPAR(XN,BN,X,B1,NPP,ITP) 37 : 38 : C-------------------------------------------------------- 39 : C DEFINITION DES DIFFERENTS DEBITS 40 : C-------------------------------------------------------- 41 : C debit air 42 : QA=QAE/E1/B1 43 : 44 : C debit eau (vapeur+liquide) 45 : QE1=QEE1/E1/B1 46 : 47 : C debit total 48 : Q=QA+QE1 49 : C 50 : C DEFINITION DES MASSES VOLUMIQUES 51 : C CONSTITUANT: 52 : C : air 53 : C : vapeur 54 : 55 : PV1=PHI1*P1 56 : 57 : C masse volumique de air 58 : ROA1=(P1-PV1)/RA/T1 59 : C WRITE(6,*) 'T1,P1,PV1,RA',T1,P1,PV1,RA 60 : 61 : C masse volumique vapeur 62 : ROV1=ROVAP0(PV1,T1) 63 : 64 : C masse volumique du melange 65 : RO1 = ROA1 + ROV1 66 : C write(6,*) ' ROA ROV RO ',ROA1,ROV1,RO1 67 : 68 : C-------------------------------------------------------- 69 : C CALCUL DE LA VITESSE DE ECOULEMENT EN ENTREE 70 : C-------------------------------------------------------- 71 : 72 : C calcul de la vitesse 73 : U1=QA/ROA1 74 : C write(6,*) ' U1 ',U1 75 : 76 : 77 : C-------------------------------------------------------- 78 : C CALCUL DU COEFFICIENT DE COMPRESSIBILITE Z de vapeur 79 : C CONSTITUANT: 80 : C : vapeur 81 : C-------------------------------------------------------- 82 : C calcul du ceofficient de compressibilite Z vapeur 83 : ZV = ZVAP0(ROV1,T1) 84 : C WRITE(6,*) 'Z',ZV 85 : 86 : 87 : C-------------------------------------------------------- 88 : C CALCUL DES CONDUCTIVITES THERMIQUES 89 : C CONSTITUANT: 90 : C : air 91 : C-------------------------------------------------------- 92 : 93 : C calcul la conductivite de air 94 : XLA = BLA(T1,T0) 95 : 96 : 97 : 98 : C-------------------------------------------------------- 99 : C CALCUL DES CHALEURS SPECIFIQUES 100 : C CONSTITUANT: 101 : C : vapeur 102 : C---------------------------------------------------- 103 : 104 : CV1 = DHVDT0(PV1,T1) 105 : C WRITE(6,*) 'cv1',CV1 106 : 107 : 108 : 109 : C-------------------------------------------------------- 110 : C CALCUL DES VISCOSITES DYNAMIQUES 111 : C CONSTITUANT: 112 : C : air + vapeur 113 : C-------------------------------------------------------- 114 : 115 : XMU=BMUG(T1,PHI1,T0) 116 : C WRITE(6,*) 'xmu',XMU 117 : 118 : 119 : C-------------------------------------------------------- 120 : C CALCUL DU NOMBRE DE PRANDT DE ECOULEMENT 121 : C CONSTITUANT: 122 : C : air + vapeur 123 : C-------------------------------------------------------- 124 : 125 : PR=XMU*CV1/BLA(T1,T0) 126 : C WRITE(6,*) 'PR',PR 127 : C-------------------------------------------------------- 128 : C CALCUL DES ENTHALPIES SPECIFIQUES 129 : C CONSTITUANT: 130 : C : eau liquide 131 : C : vapeur 132 : C-------------------------------------------------------- 133 : HVS01 = HVS0(PV1,T1) 134 : HLS01 = HLS0(P1,T1) 135 : C WRITE(6,*) 'HV',HVS01 136 : C WRITE(6,*) 'HL',HLS01 137 : 138 : C-------------------------------------------------------- 139 : C CALCUL CHALEUR LATENTE DU MELANGE 140 : C (utile pour la condensation en film) 141 : C CONSTITUANT: 142 : C : eau+vapeur 143 : C-------------------------------------------------------- 144 : 145 : XLAT = (HVS01-HLS01) 146 : C WRITE(6,*) 'L',XLAT 147 : 148 : C-------------------------------------------------------- 149 : C ACTUALISATION DES VARIABLES 150 : C CONSTITUANT: 151 : C : vitesse U 152 : C : pression totale P 153 : C : temperature T 154 : C---------------------------------------------------- 155 : 156 : E=E1 157 : B=B1 158 : U=U1 159 : T=T1 160 : P=P1 161 : RO=RO1 162 : 163 : 164 : IF ((X+DX).GT.1.) DX=1-X 165 : X=X+DX 166 : 167 : DU=0 168 : T2=0 169 : P2=0 170 : 170 : CALL BTPAR(XN,EN,X,E2,NPP,ITP) 170 : CALL BTPAR(XN,BN,X,B2,NPP,ITP) 171 : 172 : NITER=0 173 : 174 : 10 CONTINUE 175 : 176 : NITER=NITER+1 177 : IF(KIMP.GE.1) THEN 178 : write(*,*) 'NITER= ',NITER 179 : ENDIF 180 : 181 : TT2=T2 182 : PP2=P2 183 : 184 : C-------------------------------------------------------- 185 : C CALCUL DU NOMBRE DE REYNOLDS 186 : C CONSTITUANT: 187 : C : melange 188 : C------------------------------------------------------- 189 : RE=Q*2*E/XMU 190 : IF(KIMP.GE.1) THEN 191 : WRITE(6,*) 'bsur2: RE',RE 192 : WRITE(6,*) 'bsur2: Q,E,B,XMU',Q,E,B,XMU 193 : ENDIF 194 : C-------------------------------------------------------- 195 : C CALCUL DE LA PERTE DE CHARGE DP DUE AU FROTTEMENT 196 : C CONSTITUANT: 197 : C : Melange 198 : C-------------------------------------------------------- 199 : 200 : DPP = XL*Q*Q/4/E 201 : BK = BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3, 202 : *XKUT4) 203 : DP=DX*BK*DPP/RO 204 : C WRITE(6,*) 'XL Q E ', XL,Q,E 205 : C WRITE(6,*) 'DX RO DP', DX,RO,DP 206 : 207 : P2=P1-DP 208 : IF (P2.LT.PSLIM) THEN 209 : PSQ = -1. 210 : RETURN 211 : ENDIF 212 : 213 : 214 : C-------------------------------------------------------- 215 : C CALCUL DU TERME SOURCE S POUR LE CALCUL DE T 216 : C CONSTITUANT: 217 : C : melange 218 : C---------------------------------------------------- 219 : 220 : S=0.D0 221 : 222 : C-------------------------------------------------------- 223 : C CALCUL DU COEFFICIENT ECHANGE THERMIQUE H 224 : C CONSTITUANT: 225 : C : Melange 226 : C---------------------------------------------------- 227 : 228 : HM = BHECH(RE,REL,PR,XLA,E) 229 : C write(6,*) RE,REL,PR 230 : C write(6,*) XLA,E 231 : C write(6,*) ' HM ',HM 232 : 233 : C if vapour at the inlet 234 : 235 : IF(PHI1.GT.1.D-2) THEN 236 : C 237 : C wall condensation 238 : D=2.5D-5 239 : A=2.3D-5 240 : 241 : C-------------------------------------------------------- 242 : C CALCUL PRESSION DE SATURATION 243 : C CONSTITUANT: 244 : C : vapeur 245 : C---------------------------------------------------- 246 : 247 : PP=PSATT0(TP) 248 : C WRITE(6,*) 'TP',TP 249 : C-------------------------------------------------------- 250 : C CALCUL PRESSION PARTIELLE DE VAPEUR 251 : C CONSTITUANT: 252 : C : vapeur 253 : C---------------------------------------------------- 254 : C WRITE(6,*) 'ROV1',ROV1 255 : PV=PVAP0(ROV1,T1) 256 : C WRITE(6,*) 'PV',PV 257 : 258 : 259 : C-------------------------------------------------------- 260 : C CONDENSATION EN FILM 261 : C-------------------------------------------------------- 262 : 263 : 264 : RAPP=(P1-PP)/(P1-PV) 265 : 266 : IF (RAPP.GT.(1.0001)) THEN 267 : XJ=XW*(HM/CV1)*(D/A)**(2./3.)*DLOG(RAPP) 268 : ELSE 269 : XJ=0.D0 270 : ENDIF 271 : 272 : DQEE=-2.*XJ*DX*XL*B 273 : 274 : SW=(XLAT/Q/CV1/E)*DQEE 275 : 276 : DELTAT = T-TP 277 : IF (ABS(DELTAT).GT.1D-3) THEN 278 : HW=XJ*XLAT/(T-TP) 279 : C write(6,*) ' T TP HW ',T,TP,HW 280 : ENDIF 281 : 282 : QEE2=QEE1+DQEE 283 : QE2=QEE2/E2/B2 284 : C write(6,*) ' QE1 QE2 ',QE1,QE2 285 : 286 : IF ((QEE2/QEE1).LT.(-1.D-5)) THEN 287 : C write(6,*) ' relative water flowrate ',QE2/QE1 288 : RINDEX = 0.5 289 : RETURN 290 : C PHI2=0.D0 291 : C QEE2=0.D0 292 : C QE2=0.D0 293 : ELSE 294 : A=RA*QA/RV/(QE2*ZV) 295 : PHI2=1./(1.+A) 296 : ENDIF 297 : 298 : 299 : * no vapour at the inlet 300 : ELSE 301 : XLAT=(HVS01-HLS01) 302 : HW=0.D0 303 : PHI2=0.D0 304 : QEE2=0.D0 305 : QE2=0.D0 306 : ENDIF 307 : 308 : 309 : C write(6,*) ' QE ',QE1,QE2 310 : C-------------------------------------------------------- 311 : C CALCUL DES COEFFICIENTS ALPHAi 312 : C CONSTITUANT: 313 : C : on se sert des fonctions de VARI 314 : C-------------------------------------------------------- 315 : 316 : 317 : C definition des alphai pour resolution systeme 318 : alpha1=DHVDT0(PV1,T1) 319 : alpha2=DHVDP0(PV1,T1) 320 : 321 : alpha6=DZVDP0(PV1,T1) 322 : alpha7=DZVDT0(PV1,T1) 323 : alpha8=DRVDP0(PV1,T1) 324 : alpha9=DRVDT0(PV1,T1) 325 : 326 : C WRITE(6,*) 'ALPHA',alpha1 327 : C WRITE(6,*) 'ALPHA',alpha2 328 : C WRITE(6,*) 'ALPHA',alpha3 329 : C WRITE(6,*) 'ALPHA',alpha4 330 : C WRITE(6,*) 'ALPHA',alpha5 331 : C WRITE(6,*) 'ALPHA',alpha6 332 : C WRITE(6,*) 'ALPHA',alpha7 333 : C WRITE(6,*) 'ALPHA',alpha8 334 : C WRITE(6,*) 'ALPHA',alpha9 335 : 336 : C-------------------------------------------------------- 337 : C DEFINITION DU Qxx,Qy 338 : C CONSTITUANT: 339 : C : on se sert des fonctions de VARI 340 : C-------------------------------------------------------- 341 : 342 : Qxx= QAE*CA+QEE1*alpha1 343 : s1=ZV*((alpha9*T1)+ROV1*T1*alpha7) 344 : s2= RA*T1*(alpha8*ZV+ROV1*alpha6) 345 : s3= 1-s2 346 : s4=s1/s3 347 : Qy= Qxx + QEE1*alpha2*RA*s4 348 : 349 : C-------------------------------------------------------- 350 : C CALCUL DE LA TEMPERATURE T2 EN SORTIE DU TRONCON 351 : C---------------------------------------------------- 352 : 353 : HEFF = 2 * HM 354 : C bt3 traite l echange de chaleur par unite de largeur 355 : Qy = Qy/B 356 : T2 = BT3(T1,TP,S,HEFF,Qy,DX,XL) 357 : IF(KIMP.GE.1) THEN 358 : WRITE(6,*) 'bsur2 : T2= ',T2 359 : ENDIF 360 : 361 : C-------------------------------------------------------- 362 : C CALCUL DE LA VARIATION DE TEMPERATURE SUR TRONCON 363 : C-------------------------------------------------------- 364 : 365 : 366 : DT = T2-T1 367 : 368 : C-------------------------------------------------------- 369 : C CALCUL DE LA VITESSE DE ECOULEMENT U2 SORTIE TRONCON 370 : C--------------------------------------------------------- 371 : 372 : C WRITE(6,*) 'P2 RA T2 ',P2,RA,T2 373 : 374 : C initialisation avec la valeur a l'entree du troncon 375 : 376 : ZTRAN = ZV 377 : NZ = 0 378 : 20 CONTINUE 379 : NZ = NZ + 1 380 : 381 : AA = (QEE2/QAE)* (RV*ZTRAN/RA) 382 : PHI2 = AA/(1.+AA) 383 : PV2 = PHI2 * P2 384 : ROV2 = ROVAP0(PV2,T2) 385 : Z2 = ZVAP0(ROV2,T2) 386 : C write(6,*) ' AA PHI2 ROV2 ',AA,PHI2,ROV2 387 : DELTAZ = ABS(Z2 - ZTRAN) 388 : C WRITE(6,*) ' NZ ZVAP ', NZ,Z2 389 : ZTRAN=Z2 390 : IF ((DELTAZ.GT.1D-5).AND.(NZ.LE.20)) GOTO 20 391 : 392 : IF (KIMP.GE.2) WRITE(6,*) ' NZ ZVAP ', NZ,Z2 393 : 394 : ROA2 = (1.-PHI2)*P2/RA/T2 395 : RO2 = ROA2 + ROV2 396 : 397 : U2 = QA/ROA2 398 : 399 : C WRITE(6,*) ' ROA ROV RO Z PV ',ROA2,ROV2,RO2,Z2,PV2 400 : C WRITE(6,*) ' PHI U ',PHI2,U2 401 : 402 : 403 : IF(KIMP.GE.2) THEN 404 : write(6,1200) X,(P2/P0),(T2-T0),PV2,Z2 405 : 1200 FORMAT(1X,' bsur2 X P2 T2 PV2 Z2 ',5E12.5) 406 : write(6,1300) QE1,QE2 407 : 1300 FORMAT(1X,' QE1 QE2 ',2E12.5) 408 : ENDIF 409 : 410 : 411 : C--------------------------------------------------------- 412 : C CALCUL DES VALEURS MOYENNES DE P,T,U,RO 413 : C--------------------------------------------------------- 414 : 415 : T=(T1+T2)/2 416 : P=(P1+P2)/2 417 : RO= (RO1+RO2)/2 418 : C write(6,*) ' RO1 RO2 ',RO1,RO2 419 : U=Q/RO 420 : 421 : C write(6,*) ' T P RO U ',T,P,RO,U 422 : 423 : QE = (QE1 + QE2)/2 424 : Q = QA + QE 425 : 426 : QW2=(QEE2-QEE1)/DX/XL/B 427 : 428 : ERT=DABS((T2-TT2)/T2) 429 : ERP=DABS((P2-PP2)/P2) 430 : 431 : IF(((ERT.GT.1E-4).OR.(ERP.GT.1E-4)) 432 : $ .AND.(NITER.LE.10)) GOTO 10 433 : 434 : 435 : IF(KIMP.GE.2) THEN 436 : write(6,*) ' HM HW HW/HM ',HM,HW,HW/HM 437 : ENDIF 438 : 439 : C calcul des puissances echangees 440 : 441 : C RAPH=(CV1*T+XLAT0)/XLAT 442 : C H=HM+HW*RAPH 443 : H = HM 444 : 445 : 446 : C------------------------------------------- 447 : C BILAN ENERGIE TOTALE 448 : C------------------------------------------- 449 : 450 : HH1=(QAE*CA*T1) + (QEE1*HVS0(PV1,T1)) 451 : HH2=(QAE*CA*T2) + (QEE2*HVS0(PV2,T2)) 452 : C puissance fluide 453 : C ! PF = HH2-HH1 454 : PF = (HH2-HH1) - (DQEE*HVS0(PV1,T1)) 455 : 456 : C puissance paroi 457 : PP = H*(T-TP)*2*DX*XL*B 458 : 459 : 460 : IF(KIMP.GE.2) THEN 461 : write(6,1000) X,HH1,HH2,PF,PP 462 : 1000 FORMAT(1X,' bsur2 X HH1 HH2 DH HDT ',5E12.5) 463 : ENDIF 464 : IF(KIMP.NE.0.AND.X.EQ.1.) THEN 465 : write(6,2110) X,RE,BK 466 : 2110 format(1X,'bsur2 X RE BKRO ',3E12.5) 467 : ENDIF 468 : DPF=PF/DX/XL/B 469 : DPP=PP/DX/XL/B 470 : 471 : RETURN 472 : END 473 : 474 : 475 : 476 : 477 : 478 :
© Cast3M 2003 - All rights reserved.
Disclaimer